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ABSTRACT 

We present the evolution of the stellar mass function (SMF) of galaxies from z = 4.0 to z = 1.3 measured 
from a sample constructed from the deep NIR MUS YC, the FIRES, and the GOODS-CDFS surveys, all having 
very high-quality optical to mid-infrared data. This sample, unique in that it combines data from surveys with a 
large range of depths and areas in a self-consistent way, allowed us to 1) minimize the uncertainty due to cosmic 
variance and empirically quantify its contribution to the total error budget; 2) simultaneously probe the high- 
mass end and the low-mass end (down to ^ 0.05 times the characteristic stellar mass) of the SMF with good 
statistics; and 3) empirically derive the redshift-dependent completeness limits in stellar mass. We provide, for 
the first time, a comprehensive analysis of random and systematic uncertainties affecting the derived SMFs, 
including the effect of metallicity, extinction law, stellar population synthesis model, and initial mass function. 
We find that the mass density evolves by a factor of ^ 17^^, since z = 4.0, mostly driven by a change in the 
normalization $*. If only random errors are taken into account, we find evidence for mass-dependent evolution, 
with the low-mass end evolving more rapidly than the high-mass end. However, we show that this result is no 
longer robust when systematic uncertainties due to the SED-modeling assumptions are taken into account. 
Another significant uncertainty is the contribution to the overall stellar mass density of galaxies below our 
mass limit; future studies with WFC3 will provide better constraints on the SMF at masses below lO"* Mq at 
z > 2. Taking our results at face value, we find that they are in conflict with semi-analytic models of galaxy 
formation. The models predict SMFs that are in general too steep, with too many low-mass galaxies and too 
few high-mass galaxies. The discrepancy at the high-mass end is susceptible to uncertainties in the models and 
the data, but the discrepancy at the low-mass end may be more difficult to explain. 

Subject headings: galaxies: distances and redshifts — galaxies: evolution — galaxies: formation — galaxies: 
fundamental parameters — galaxies: high-redshift — galaxies: luminosity function, mass 
function — galaxies: stellar content — infrared: galaxies 



1. INTRODUCTION 

Understanding the formation mechanisms and evolution 
with cosmic time of galaxies is one of the major goals of ob- 
servational cosmology. In the current picture of structure for- 
mation, dark matter halos build up in a hierarchical fashion 
through the dissipationless mechanism of gravitational insta- 
bility controlled by the nature of the dark matter, the power 
spectrum of density fluctuations, and the parameters of the 
cosmological model. The assembly of the stellar content of 
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galaxies is governed by much more complicated physics, such 
as the mechanisms of star formation, gaseous dissipation, the 
feedback of stellar and central super-massive black hole en- 
ergetic output on the baryonic material of the galaxies, and 
mergers (see Baugh 2006, and references therein for a primer 
on hierarchical galaxy formation). 

A powerful approach to understand these physical pro- 
cesses (including their relative importance as function of cos- 
mic time) is to directly witness the growth of the stellar con- 
tent in galaxies. Galaxies can grow their stellar mass both 
from in-situ star formation and/or merger events. Determin- 
ing the growth of their stellar content as a function of both 
redshift and stellar mass provides insights into the physical 
processes governing the assembly and the evolution of galax- 
ies. The stellar mass function (SMF) of galaxies and its evo- 
lution with cosmic time represents therefore a powerful tool 
to directly measure the build-up of the stellar mass content of 
galaxies. 

In the past decade, significant observational progress has 
been made in the measurement of the SMF of galaxies and 
its evolution with redshift. Locally, the SMF has been mea- 
sured from the 2dF Galaxy Redshift Survey (Cole et al. 2001) 
and the Sloan Digital Sky Survey (SDSS; Bell et al. 2003), 
providing the z ~ benchmark. At intermediate redshifts 
(z < 1 .4), the SMF has also been measured to a satisfactory 
degree from the VIMOS VLT Deep Survey (VVDS; Pozzetti 
et al. 2007; Vergani et al. 2008), the DEEP-2 Galaxy Red- 
shift Survey (Bundy et al. 2006), and the COMBO-17 survey 
(Borch et al. 2006). The SMF appears to evolve slowly at 
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z < 1, with about half of the total stellar mass density at z ~ 
already in place at z ^ 1 . 

The SMF has also been measured at higher redshifts, up 
to z ^ 5, using photometric redshifts derived from multi- 
waveband imaging surveys (e.g., Dickinson et al. 2003; 
Fontana et al. 2003; Drory et al. 2004; Fontana et al. 2004; 
Conselice et al. 2005 Drory et al. 2005; Fontana et al. 2006; 
Eisner et al. 2008; Perez-Gonzalez et al. 2008). The general 
consensus is that at z > 1 the stellar mass assembly proceeds 
much more quickly than at lower redshifts. The density of 
massive galaxies appears to have strongly evolved in the red- 
shift range 1 .5 < z <^ 3 (e.g., Fontana et al. 2006), with roughly 
40% of the local stellar mass density assembled between z = 4 
andz = 1 (e.g., Drory et al. 2005; Perez-Gonzalez et al. 2008). 

However, at redshift z > 2, the derived SMFs are not in very 
good agreement (e.g., Eisner et al. 2008), possibly caused by 
differences in modeling techniques, field-to-field variations, 
and/or systematic differences in photometric redshift estimate 
errors. The cause of these differences is difficult to isolate, 
partly because most existing studies are characterized by ei- 
ther wide but very shallow data (e.g., Drory et al. 2004), or 
by deep data over relatively small areas (e.g., Dickinson et 
al. 2003; Fontana et al. 2003; Consehce et al. 2005; Drory et 
al. 2005; Fontana et al. 2006; Eisner et al. 2008). Therefore, 
they either only probe the high-mass end of the high-z SMF, 
or they survey volumes of the universe small enough that cos- 
mic variance constitutes a significant, if not dominant, source 
of uncertainties. Only the IRAC-selected sample constructed 
by Perez-Gonzalez et al. (2008) provides relatively deep data 
over three fields for a total of ~660 arcmin^. Furthermore, 
a comprehensive analysis of the errors affecting the derived 
SMFs at high redshift is still missing in the literature. Often, 
only Poisson errors are considered. Uncertainties from photo- 
metric redshift errors are seldomly included, and cosmic vari- 
ance is never included (which, as we will show, is a very large, 
if not dominant, component of random errors even with rel- 
atively large surveyed area of ^600 arcmin^). Perhaps even 
more important is the missing analysis of systematic uncer- 
tainties in the derived SMF caused by the assumptions in the 
SED modehng to estimate the stellar masses. For example, 
a satisfactory agreement between different stellar population 
synthesis models has yet to be achieved, and there is evidence 
for evolution with redshift of some of the relevant parameters, 
such as the IMF (e.g., Dave 2008; van Dokkum 2008; Wilkins 
et al. 2008). 

In this paper we take advantage of the very high-quality 
data from the optical to the mid-infrared (MIR) available over 
the deep NIR MUSYC, the GOODS-CDFS, and the FIRES 

surveys to measure the evolution of the SMF of galaxies from 
z = 4.0 to z = 1.3 and to provide the first comprehensive anal- 
ysis of its random and systematics uncertainties. The com- 
posite A'-selected sample, constructed from a total surveyed 
area of ~600 arcmin-^ and probing as deep as K^'\5 a) = 25.6 
, is unique in that it combines data from surveys with a large 
range of depths and areas in a self-consistent way. Crucially, 
the data handhng, creation of catalogs, and determination 
of photometric redshifts and masses were all done with the 
same methods and techniques. The composite sample there- 
fore self-consistently combines the advantages of deep, pencil 
beam surveys with those of shallow, wide surveys, allowing us 
to: 1) minimize the uncertainties due to cosmic variance, and 
empirically quantify its contribution to the total error budget 
by exploiting the large number of independent field of views; 
2) simultaneously probe the high-mass end and the low-mass 



end of the SMF with good statistics; and 3) empirically derive 
the redshift-dependent completeness hmits by exploiting the 
different depths of the three surveys. 

This paper is structured as follows. In § 2 we present the 
composite A'-selected sample used to measure the SMF of 
galaxies at 1.3 < z < 4.0; in § 3 we describe the approach 
adopted to estimate stellar masses and the default set of SED- 
modeling assumptions. The methods used to derive the SMF 
(the l/Vmax and the maximum-likelihood methods) and the 
approach used to estimate the completeness hmits in stellar 
mass are presented in § 4, as well as the SMFs of galaxies at 
1.3 < z < 2.0, 2.0 < z < 3.0, and 3.0 < z < 4.0. A detailed 
and comprehensive analysis of the random and systematic un- 
certainties of the derived SMFs is presented in § 5, while the 
evolution of the stellar mass densities is presented in § 6. A 
comparison with the predictions from the latest generation 
of galaxy formation models is presented in § 7. Our results 
are summarized in § 8. We assume = 0.3, Qa = 0.7, and 
Ho = 70 km s"' Mpc"' throughout the paper. AH magnitudes 
are on the AB system. 

2. THE «:-SELECTED SAMPLE 

2.1. Data 

The data set we have used to estimate the SMF consists of 
a composite /iT-selected sample of galaxies built from three 
deep multi-wavelength surveys, all having very high-quality 
optical to mid-IR photometry: the "ultra-deep" Faint InfraRed 
Extragalactic Survey (FIRES; Franx et al. 2003), the Great 
Observatories Origins Deep Survey (GOODS; Giavalisco et 
al. 2004; Chandra Deep Field South [CDFS]), and the MUlti- 
wavelength Survey by Yale-Chile (MUSYC; Gawiser et al. 
2006). Photometric catalogs were created for all fields in the 
same way, following the procedures of Labbe et al. (2003). 

2.1.1. MUSYC 

The deep NIR MUSYC survey consists of four ~ 10' x 10' 
fields, namely, Hubble Deep Field-South 1 and 2 (HDFS-1, 
HDFS-2, hereafter), the SDSS-1030 field, and the CW-1255 
field, observed with the ISPl camera at the Cerro Tololo Inter- 
American Observatory (CTIO) Blanco 4 m telescope, for a to- 
tal surveyed area of ~ 430 arcmin^ (~ 416 arcmin^ not over- 
lapping). A complete description of the deep NIR MUSYC 
observations, reduction procedures, and the construction of 
the /T-selected catalog with U -to-K photometry is presented 
in Quadri et al. (2007). 

We added deep Spitzer-IRAC 3.6-8.0 /um imaging to the 
NIR MUSYC survey. The IRAC data over the HDFS-1 are 
part of the GTO-214 program (PL: Fazio), while the IRAC 
data over the other three fields come from the GO-30873 pro- 
gram (PL: Labbe). The average total limiting magnitudes 
of the IRAC images are ~24.5, 24.2, 22.4, and 22.3 (3 a, 
AB magnitude) in the 3.6, 4.5, 5.8, and 8.0 /xm bands, re- 
spectively.^ The IRAC data, reduction, and photometry are 
described in detail in Appendix A. 

The /T-selected catalogs with IRAC photometry included 
is publicly available at http://www.astro.yale.edu/musyc. The 
SDSS-1030, CW-1255, HDFS-1, and HDFS-2 catalogs are Kg 
band-limited multicolor source catalogs down to K^^ =23.6, 
23.4, 23.7, and 23.2, for a total of 3273, 2445, 2996, and 21 18 

' For tlie MUSYC survey, less tlian 4% of ttie A'-selected sources tiave 
IRAC S/N < 3 in tlie 3.6 band; for tlie FIRES and FIREWORKS sam- 
ples, which are significantly deeper than the MUSYC sample, ~8% of the 
galaxies have IRAC S/N < 3. 
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sources, over fields of ~ 109, ~ 105, ~ 109, ~ 106 arcmin^, 
respectively. All four fields were exposed in 14 different 
bands, U, B, V, R, I, z, J, H, K, and the four IRAC chan- 
nels. The SDSS-1030, CW-1255, HDFS-1, and HDFS-2 K- 
selected catalogs have 90% completeness levels at K^^ =23.2, 
22.8, 23.0 and 22.7, respectively. The final catalogs used in 
the construction of the composite sample have 2825, 2197, 
2266, and 1749 objects brighter than the 90% completeness 
in the Kg band, over an effective area of 98.2, 91.0, 97.6, and 
85.9 arcmin^, respectively, for a total of 9037 sources over 
372.7 arcniin^. 

2.1.2. FIRES 

FIRES consists of two fields, namely, the Hubble Deep 
Field-South proper (HDF-S) and the field around MS 1054- 
03, a foreground cluster at z = 0.83. A complete descrip- 
tion of the FIRES observations, reduction procedures, and 
the construction of photometric catalogs is presented in detail 
in Labbe et al. (2003) and Forster Schreiber et al. (2006) for 
HDF-S and MS 1054-03 (hereafter HDFS and MS-1054, re- 
spectively). Both /Ts-selected catalogs were later augmented 
with Spitzer-IRAC data (Wuyts et al. 2007; Toft et al. 2007). 
The HDFS catalog has 833 sources down to K^'^ =26.0 over 
an area of 2.5' x 2.5'. The MS-1054 catalog has 1858 sources 
down to K'°^ =25.0 over an area of 5.5' x 5.3'. The HDFS field 
was exposed in the WFPC2 U^qo, B450, Veoe, hi4 pass-bands, 
the ISAAC Js, H, and Ks bands, and the four IRAC chan- 
nels. The MS-1054 /Ts-selected catalog comprises FORSl U, 
B, V, WFPC2 V(m and hu, ISAAC /, H, and K^, and IRAC 
3.6 //,m - 8.0 /im photometry. The HDFS and MS-1054 cat- 
alogs have 90% completeness levels al /Tg"' =25.5 and 24.1, 
respectively. The final HDFS and MS-1054 catalogs used in 
the construction of the composite sample have 715 and 1547 
objects brighter than the 90% completeness in the band, 
over an effective area of 4.5 and 21.0 arcmin-^, respectively. 

2.1.3. FIREWORKS-CDFS 

In this work, we adopted the /sTs-selected catalog (dubbed 
FIREWORKS) of the CDFS field consfincted based on the 
pubhcly available GOODS-CDFS data by Wuyts et al. (2008). 
The photometry was performed in an identical way to that 
of the FIRES fields, and the included passbands are the ACS 
B435, V'eoe, ins, and zgso bands, the WFI f/ag, B, V, R, and / 
bands, the ISAAC J, H, and Ks bands, and the 4 IRAC chan- 
nels. The /iTs-selected catalog comprises 6308 objects down 
to 7^5°' =24.6 over a total surveyed area of 138 arcmin^; the 
variation in exposure time and observing conditions between 
the different ISAAC pointings lead to an in-homogeneous 
depth over the whole GOODS-CDFS field (hereafter, CDFS). 
The final CDFS catalog used in the construction of the com- 
posite sample comprises 3559 objects brighter than the 90% 
completeness level (K^* = 23.7), over an effective area of 
113 arcmin^ with coverage in all bands. 

2.2. Photometric redshifts 

The fraction of isT-selected galaxies with spectroscopic red- 
shifts is ~10% (6%) with Zspec > (Zspec ^ 1-3) in the com- 
posite /T-selected sample. Consequently, we must rely pri- 
marily on photometric redshift estimates. Photometric red- 
shifts Zphot for all galaxies were derived using the EAZY pho- 
tometric redshift code (Brammer et al. 2008). EAZY fits the 
observed SED of each galaxy with a non-negative linear com- 
bination of galaxy templates. The template set used in this 
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Fig. 1. — Spectroscopic vs. photometric redsliifts for tlie deep NIR 
MUSYC (open triangles), the CDFS FIREWORKS {gray filled circles), 
the MS-1054 {open squares) and the HDFS {open stars) FIRES surveys, 
shown on a "pseudo-log" scale. The total number N of Zspec, the median in 
A2/(l +zspec), the (Tnmad. and the number of catastrophic outliers is speci- 
fied (values in parenthesis refer to galaxies with Zspec ^1.3). The comparison 
between Zspcc and Zphot is extremely good, both at low and at high redshifts. 

work is the default EAZY template set eazy_vl.O. '"^ This 
template set was constructed from a large number of PEGASE 
models (Fioc & Rocca-Volmerange 1997). It consists of 5 
principal component templates that span the colors of galax- 
ies in the semi-analytic model by De Lucia & Blaizot (2007), 
plus an additional template representing a young (50 Myr) 
and heavily obscured (Av=2.75) stellar population to account 
for the existence of dustier galaxies than present in the semi- 
analytic model. The default template error function (TEM- 
PLATE_ERROR.eazy_vl.O) was applied to down-weight the 
rest-frame UV and rest-frame NIR during the fitting proce- 
dure. The /T-band magnitude was used as a prior in construct- 
ing the redshift probabiUty distribution for each galaxy, and 
the value Zpm (the redshift marginalized over the total prob- 
ability distribution) was adopted as the best estimate of the 
galaxy redshift. 

Figure 1 shows the comparison of Zphot versus Zspec for 
the /T-selected samples. For the full sample, the median in 
Az/ (1 +Zspec), with Az = Zphot-Zspec is 0.000, with the normal- 
ized median absolute deviation ctnmad" =0.033. The frac- 
tion of catastrophic outliers, here defined as galaxies with 
Az/(1 +Zspec) > 5ctnmad, is 4%. For galaxies with Zspec > 1 -3, 
the median in Az/(1 +Zspec) is -0.015, with ctnmad = 0.061, 
and the fraction of catastrophic outliers is 5% (9% if (tnmad = 
0.033 corresponding to the entire sample is used in the defini- 
tion of catastrophic outhers). 

Restricting the comparison between Zphot and Zspec for the 

The template set needs to be large enough that it spans the broad range 
of multi-band galaxy colors and small enough that the color and redshift de- 
generacies are kept to a minimum (e.g., Benitez 2000). The default template 
set used in this work to estimate photometric redshifts was carefully con- 
structed and tested in Brammer et al. (2008). It has been shown to satisfy the 
requirements for a satisfactory tempkite sets, providing significantly reduced 
systematic effects and smaller scatter in the z^hot vs Zspcc at all redshifts. 

" The normalized median absolute deviation ctnmad. defined as 1.48 X 
median[\{J\7.-tnedian{/\z))l(i +Zspec)|], is equal to the standard deviation 
for a Gaussian distribution, and it is less sensitive to outliers than the usual 
definition of the standard deviation (e.g. Ilbert et al. 2006) 
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MUSYC sample alone gives similarly good results. We note 
that the estimate of the photometric redshifts in MUSYC has 
improved with respect to previously published works, due to 
the use of EAZY and the inclusion of the IRAC data. 

The effects of photometric redshift errors (both random and 
systematic) on the derived SMFs and stellar mass densities are 
quantified and discussed in § 5.1. Specifically, the sensitivity 
of the derived SMFs to the choice of template set is tested in 
§5.1.2, and its contribution included in the total error budget. 

2.3. The composite K-selected sample 

Stars in all /T-selected catalogs were identified by spec- 
troscopy, by fitting the object SEDs with stellar templates 
from Hauschildt, Allard, & Baron (1999) and/or inspecting 
their morphologies, as in Rudnick et al. (2003). On average, 
approximately 12% of all objects were classified as stars. 

We constructed a composite /T-selected sample of high- 
redshift (1.3 < z < 4.0) galaxies to be used in the derivation 
of SMFs of galaxies in § 4. The large surveyed area of the 
deep NIR MUSYC and FIREWORKS -CDFS surveys allows 
us to probe the high-mass end of the SMF with good statis- 
tics, as well as empirically quantify uncertainties due to field- 
to-field variations. The very deep FIRES survey is critical in 
allowing us to constrain the low-mass end of the SMF. Finally, 
the FIREWORKS -CDFS survey bridges the two sUghtly over- 
lapping regimes probed by the MUSYC and the FIRES sur- 
veys. The final composite sample includes 2014, 830, and 
213 A'-selected galaxies in the three targeted redshift inter- 
vals 1.3 < z < 2.0, 2.0 < z < 3.0, and 3.0 < z < 4.0, respec- 
tively, for a total of 3057 galaxies over an effective area of 
511.2 arcmin^ with f^" < 25.5 at 1.3 < z < 4.0. Of these, 
~6% have spectroscopic redshifts. 

3. SED MODELING: DEFAULT ASSUMPTIONS 

Stellar masses for the HDFS, the MS-1054, and the CDFS 
fields were derived in Forster Schreiber et al. (in prep.), 
following the procedure described in Wuyts et al. (2007). 
We adopted the same procedure to derive the stellar masses 
of the galaxies in the MUSYC HDFS-1, HDFS-2, SDSS- 
1030, and CW-1255 fields. In the following we describe our 
default assumptions, {BC03,ZQ,Kroupa,Calzetti), to perform 
SED modeling. 

We have generated stellar population synthesis models with 
the evolutionary synthesis code developed by Bruzual & 
Chai-lot (2003) (BC03). We selected the "Padova 1994" evo- 
lutionary tracks (Fagotto et al. 1994), which are preferred 
by Bruzual & Chariot over the more recent "Padova 2000" 
tracks because the latter may be less reliable and predict a 
hotter red giant branch leading to worse agreement with ob- 
served galaxy colors. The solar metallicity set of tracks was 
used. We fitted the BC03 templates to the observed optical-to- 
8 yum SED with the HYPERZ stellar population fitting code, 
version 1.1 (Bolzonella et al. 2000). We allowed the fol- 
lowing star formation histories: a single stellar population 
(SSP) without dust, a constant star formation history (CSF) 
with dust, and an exponentially declining star formation his- 
tory with an e-folding timescale of 300 Myr (T300) with dust. 
The allowed Ay values ranged from to 4 in step of 0.2 mag, 
and we used the attenuation law of Calzetti et al. (2000), de- 
rived empirically from observations of local UV-bright star- 
burst galaxies under the formalism of a foreground screen of 
obscuring dust. We constrained the time since the onset of star 
formation to lie between 50 Myr and the age of the universe 
at the respective redshift. Finally, we scaled from a Salpeter 
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Fig. 2. — Stellar mass versus redshift for FIRES HDFS (red filled circles), 
FIRES MS-1054 (blue open squares), FIREWORKS CDFS (orange open tri- 
angles), and MUSYC (green filled circles). 

(1955) IMF with lower and upper mass cutoffs of 0.1 Mq 
and 100 Mq to a pseudo-Ki'oupa (2001) IMF by dividing the 
stellar masses by a factor of 1.6. Therefore, the adopted de- 
fault SED-modeling assumptions can be summarized with the 
combination {BC03,ZQ,Kroupa, Calzetti), where the first, sec- 
ond, third, and fourth elements are the adopted stellar popu- 
lation synthesis model, the metallicity, the IMF, and the ex- 
tinction law, respectively. Note that we have used the mass 
of living stars plus stellar remnants, Mstar, instead of the total 
mass of stars formed, Mtot. Mtot is the integral of the SFR and 
corresponds to the total gas mass consumed since the onset of 
star formation, whereas Mstar is computed by subtracting from 
Mtot the mass returned to the ISM by evolved stars via stellar 
winds and supernova explosions (for details see Bruzual & 
Chai-lot 2003). 

In Figure 2 we show the stellar mass Mstar versus the red- 
shift for the /T-selected composite sample in the targeted red- 
shift range 1.3 < z < 4.0. 

Whereas (BC03,ZQ,Kroupa,Calzetti) represents our default 
SED-modeling assumptions, we have considered different 
stellar population synthesis models, metallicities, IMFs, and 
extinction curves to derive stellar masses. These perturbations 
on the default set of SED-modeling assumptions are described 
in detail in § 5.3, together with the resulting systematics ef- 
fects on the derived SMFs. This has generally not been ad- 
dressed in previous studies. 

One possible limitation of our approach to derive stellar 
masses in our sample is contamination by obscured AGNs, 
as they could contribute to the emission in the rest-frame NIR 
(rest-frame wavelengths longer than 2 /xm). Several recent 
studies (e.g., Papovich et al. 2006; Kriek et al. 2007; Daddi et 
al. 2007) have suggested that the fraction of obscured AGNs 
increases with redshift and stellar mass. According to Kriek et 
al. (2007), the fraction is about ^^20% for massive galaxies at 
2 < z < 2.7. To estimate the robustness of our estimated stel- 
lar masses in potential obscured AGNs, we have re-estimated 
the stellar masses of the galaxies with 8 /im excess (with 
respect to the best fit stellar population model) without in- 
cluding the two reddest IRAC channels in the SED-modeling 
(which probe the rest-frame NIR at z > 1 .3). The median dif- 
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ference of the stellar masses estimated with and without the 
5.8 and 8.0 /im IRAC bands for 8 /im excess sources is ^0%, 
while the mean difference is ^ 6%, negligible with respect 
to the typical random error (>30%; e.g., Wuyts et al. 2007; 
Muzzin et al. 2009). We therefore conclude that our stellar 
mass estimates are robust even for galaxies which likely har- 
bor an AGNs. 

4. THE STELLAR MASS FUNCTION 

In this section we describe a new technique used to properly 
quantify the completeness in stellar mass of the A'-selected 
sample used in the derivation of the SMFs (§ 4.1), and we 
discuss the two independent methods adopted to derived the 
SMFs of galaxies, namely an extended version of the 1 /Vmax 
method (§ 4.2) and the maximum- likelihood method (§ 4.3). 

The uncertainties on the derived SMFs due to photometric 
redshifts errors, cosmic variance, and different SED-modeling 
assumptions are quantified and discussed in § 5.1, § 5.2, and 
§ 5.3, respectively. 

4.1. Completeness in stellar mass 

One of the most critical steps in properly deriving the SMF 
of galaxies is to understand the completeness in stellar mass 
as function of redshift of the sample. Determining the com- 
pleteness for a flux-hmited sample is particularly challenging, 
since there is not a sharp limit in stellar mass corresponding 
to the sharp limit in flux. This is a direct consequence of the 
fact that, at any given luminosity, galaxies exhibit a range in 
mass-to-light ratios (M/L). Because of this, galaxies that are 
not selected because their observed flux is below the detec- 
tion limit can still have stellar masses well within the range 
of interest if their M/L ratios are large enough. For the same 
reason, galaxies with small stellar masses are observed in even 
relatively shallow surveys if they have small M/L ratios. 

This is obvious from the left panel of Figure 3, where 
the stellar mass Mstar is plotted versus the observed A'-band 
magnitude for the /T-selected galaxies in the narrow redshift 
range 1 .5 < z < 1 .7 (because of the narrow redshift range, 
the observed /T-band magnitude is approximately equivalent 
to the rest- frame luminosity). At a given luminosity, there is 
a range of stellar masses. The observed range in M/L ratios 
at 1.5 < z < 1.7 is largest for galaxies with intermediate lu- 
minosities and narrower for bright and for faint galaxies. For 
bright galaxies, the range in M/L ratio is smaller due to the 
fact that the bright end of the luminosity function is domi- 
nated by red galaxies with generally high M/L ratios. At the 
faint end, most galaxies are blue (e.g., Giallongo et al. 2005; 
Zucca et al. 2006; Marchesini et al. 2007) and have low M/L 
ratios, resulting in a narrow range in M/L ratios. The average 
M/L ratio for galaxies in the redshift range 1.5 < z < 1.7 at 
the faint end is a factor of ~ 5 smaller than those at the bright 
end. 

The distribution of galaxies in the left panel of Figure 3 is a 
function of redshift. First, as higher redshifts are probed, faint 
galaxies are progressively missed due to the limiting flux of 
the sample. Consequently, at higher redshifts the sample be- 
comes progressively dominated by galaxies with larger M/L 
ratios. Second, the galaxy properties will evolve with redshift, 
as well as their relative contributions in different luminosity 
ranges. For example, the bright end of the luminosity function 
is shown to be dominated by red galaxies at z ~ 2.2, whereas 
at z ^ 3 the contribution of red and blue galaxies to the ob- 
served number densities appears to be similar (Marchesini et 
al. 2007). Consequently, the observed range of M/L ratios at 



the bright end would increase going to higher redshifts. 

Most of the previous works deriving the SMFs of galaxies 
(e.g., Drory et al. 2005; Fontana et al. 2006; Perez-Gonzalez 
et al. 2008; see also Appendix C for a detailed discussion of 
previously pubUshed works and their assumed completeness 
hmits in stellar mass) have estimated their completeness in 
stellar mass based on the maximal stellar mass allowed for a 
galaxy at the flux limit of the sample. This maximal mass is 
typically taken to be the stellar mass of a passively evolving 
population with no dust extinction formed in a single burst at 
very high (z ~ 10) redshift, scaled to match the limiting flux 
(SSP-derived completeness limit, hereafter). This approach, 
although easy to implement, is conceptually wrong, and po- 
tentially affected by several caveats. 

For example, the faint-end of the galaxy luminosity func- 
tion is dominated by intrinsically blue galaxies (e.g., Zucca 
et al. 2006, Marchesini et al. 2007), characterized by small 
M/L ratios (see also left panel of Figure 3). If an SSP-derived 
completeness limit (shown by the dashed curves in the right 
panel of Figure 3) were to be used in deriving the complete- 
ness in stellar mass of a deep survey, the sample would be cut 
at conservatively too high stellar masses. More importantly, 
the derived number density at the low-mass end would be sig- 
nificantly affected due to miscalculation of Vmax in the 1 /V,nax 
method. Another compUcation is caused by the ubiquitous 
presence of dust in real galaxies. Dust extinction can reduce 
the completeness derived from no-extinction SSPs since, to 
first order, it moves the dashed curves in the right panel of 
Figure 3 upwards. While at low redshift massive galaxies 
are usually characterized by passively evolving populations 
with no or Uttle dust extinction, this does not seem to be the 
case at high redshifts, where significant amount of extinction 
(Ay ~ 1 - 2) can be present even in massive galaxies (e.g., 
Blain et al. 1999b; Muzzin et al. 2009). 

In order to avoid these problems, we have used a differ- 
ent approach to estimate the redshift-dependent completeness 
limit in stellar mass. Our approach exploits the availability 
of several samples with different depths, and the complete- 
ness of a sample is estimated empirically from the available 
deeper samples. To estimate the redshift-dependent stellar 
mass completeness Umit of one of the considered samples, 
we first selected galaxies belonging to the available deeper 
samples. Then, we scaled their fluxes and stellar masses to 
match the K-hand completeness limit of the sample we want 
to derive the completeness hmit in stellar mass for. The upper 
envelope of points in the (Mgtar.scaied-z) space, encompassing 
95% of the points, represents the most massive galaxies at the 
considered flux limit, and so provides a redshift-dependent 
stellar mass completeness limit for the considered sample. 
This method is illustrated in detail for the SDSS-1030 sam- 
ple in Appendix B. By repeating this procedure, we derived 
the redshift-dependent completeness limits in stellar mass for 
all samples. These limits are shown in the right panel of Fig- 
ure 3, along with SSP-derived completeness limits. 

Interestingly, the empirically-derived completeness limit is 
similar to the SSP-derived completeness only for the CDFS 
sample. For shallower samples (e.g., SDSS-1030), the sam- 
ple is actually less complete than what would be estimated 
with the SSP-derived completeness over the redshift range 
2.5 < z < 4. Conversely, the SSP-derived completeness is too 
conservative (by ~0.3 dex) for deeper samples, such as the 
FIRES MS- 1054 sample. 

For the MS-1054 sample, the empirically-derived com- 
pleteness at 3.0 < z < 4.0 is poorly derived due to the smaU 



Marchesini et al. 



cn 
o 




£ 10 : 



2 10 



CP 

o 



26 



24 22 
[AB] 



20 




Fig. 3. — Left panel: Stellar mass Mstai versus observed total iT-band magnitude for the galaxies at redshift 1.5 < z < 1.7 in the A'-selected sample. The thick 
black line represents the stellai' mass of a single stellar population formed at z ~ 10 and with no dust as a function of the observed S'-band magnitude; the gray 
shaded region around it represents the scatter due to the width of the considered redshift interval. The gray curves represent M/L ratios ~2, 4, 8, 16, and 32 
times smaller than that of the single stellar population formed at z ~ 10. The vertical dashed lines represent the K-hand completeness limits, from left to right, of 
the HDFS, MS-1054, CDFS, SDSS-1030, CW-1255, HDFS-1, and HDFS-2 samples. The largest range in M/L ratios is observed at intermediate luminosities, 
while the range in M/L is much narrower both at bright and faint luminosities. Right panel: Empirically-derived completeness Hmits in stellar mass as a function 
of redshift for our A'-selected samples (the HDFS-1, HDFS-2, and CW-1255 samples are omitted for clarity), plotted as solid black curves (from top to bottom, 
SDSS-1030, CDFS, MS-1054, and HDFS). For comparison, the dashed curves represent the SSP-derived completeness limits. The gray solid lines represent the 
adopted conservative completeness limit at 3.0 < z < 4.0 for the MS-1054 sample, and at all redshifts for the HDFS sample to ensure robust derivation of the 
SMFs. Note the differences, as function of redshift and stellar mass, between the empirically- and SSP-derived completeness limits. 



number of galaxies in the HDFS sample in this redshift range. 
Therefore, we have conservatively assumed its largest value 
throughout the entire redshift range 3.0 < z < 4.0. For the 
HDFS sample, there is no deeper survey that can be used 
to empirically derive its completeness in stellar mass. The 
completeness of HDFS was therefore estimated by scaling the 
empirically-derived completeness limit of the MS-1054 sam- 
ple to match the /T-band 90% completeness flux limit of the 
HDFS sample, and by taking the largest value of the scaled 
completeness within each individual redshift interval (shown 
as gray solid lines in the right panel of Figure 3). Although 
this is a very conservative approach, it ensures correct deter- 
mination of the SMF at the low-mass end. 

4.2. The 1 /Vmax method 

To estimate the observed SMF for our composite sample, 
we have applied an extended version of the 1 /Vmax algorithm 
(Schmidt 1968) as defined in Avni & Bahcafl (1980) so that 
several samples with different depths can be combined in one 
calculation. This method is described in detail in Marchesini 
et al. (2007), where it was used to derive the rest-frame opti- 
cal luminosity functions of galaxies at redshift 2.0 < z < 3.5 
from a similar /T-selected sample. The empirically-derived 
redshift-dependent completeness limits in stellar mass as de- 
rived in § 4. 1 for the individual /T-selected samples was used 
in the calculation of V,nax- The Poisson error in each stellar 
mass bin was computed adopting the recipe of Gehrels (1986). 

The 1 /Vmax estimator has the advantages of simplicity and 
no a priori assumption of a functional form for the stellar mass 
distribution; it also yields a fully normalized solution. How- 
ever, it can be affected by the presence of clustering in the 
sample, leading to a poor estimate of the faint-end slope of the 
SMF. Field-to-field variation represents a significant source 
of uncertainty in deep surveys, since they are characterized 



by small areas and hence small probed volumes. The con- 
tribution due to sample variance to the total error budget is 
quantified in § 5.2. 

4.3. The maximum likelihood method 

We also measured the observed SMF using the STY method 
(Sandage, Tammann & Yahil 1979), which is a parametric 
maximum-likelihood estimator The STY method has been 
shown to be unbiased with respect to density inhomogeneities 
(e.g., Efstathiou, Ellis & Peterson 1988), it has well-defined 
asymptotic error properties (e.g. Kendall & Stuart 1961) and 
does not require to select bin widths. 

We have assumed that the number density ^(M) of galaxies 
is described by a Schechter (1976) function, 

$(M) = (In 10)$* [lO(*^-*^'('+"'] X exp [- 10*^"^*'] , (1) 

where M = log (Mstar/M©), a is the low mass-end slope, M* = 
log(M*[jy./M0) is the characteristic stellar mass at which the 
SMF exhibits a rapid change in the slope, and $* is the nor- 
malization. 

The implementation of the STY method and the method of 
estimating errors is described in detail in Marchesini et al. 
(2007). The best-fit solution is obtained by maximizing the 
likelihood A with respect to the parameters a and M* . The 
value of is then obtained by imposing a normalization on 
the best-fit SMF such that the total number of observed galax- 
ies in the composite sample is reproduced. The 1 and 2 a 
errors on $* are estimated from the minimum and maximum 
values of $* allowed by the 1 and 2 a confidence contours in 
the (a-M*(^j) parameter space, respectively. 

4.4. Stellar Mass Functions 

Figure 4 (top panel) shows the SMFs of galaxies at redshift 
1.3 < z < 2.0, 2.0 < z < 3.0, and 3.0 < z < 4.0. Points with 
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Fig. 4. — Top panel: SMFs of galaxies at redshift 1.3 < z < 2.0 (blue), 
2.0 < z < 3.0 (green), and 3.0 < z < 4.0 (red). The filled symbols repre- 
sent the SMFs derived with the 1 /Vmax method, with error bars including 
only Poisson eiTors. The solid curves represent the SMFs derived with the 
maximum-likelihood analysis, with shaded regions representing the 1 and 
2 cr uncertainties. The an'ows show the best estimates of Mj^^ and the cor- 
responding 1 (T errors derived with the maximum-likelihood analysis. The 
black solid curve and points represent the local (z ~ 0.1) SMF from Cole 
et al. (2001). Bottom panel: (a — M*^.^^) parameter space derived from the 
maximum-likelihood analysis. Filled circles are the best-fit values of a and 
M*(^^, while the curves represent their 1 and 2 a contour levels; the colors are 
the same as in the top panel. The black filled square represent the redshift 
z ~ 0.1 value from Cole et al. (2001). Very little evolution of the shape of the 
SMF is observed from z = 3.0 to z = 1.3, and most of the evolution is in the 
characteristic density <I>*. The shape of the SMF at z = 3.5 is different, char- 
acterized by a much steeper low-mass end slope. The characteristic stellar 
mass M*,jn. seems to have evolved little, if any, from z = 2.5 to z ~ 0.1. 



error bars show the SMFs derived using the 1 /Vmax method. 
The solid curves show the SMFs derived with the maximum- 
likelihood analysis, while the shaded regions represent their 
1 and 2 a uncertainties. The bottom panel of Figure 4 shows 
the best-fit value and the 1 and 2 a confidence contour lev- 
els of the two Schechter function parameters a and M*^.^^ in 
the three targeted redshift intervals. The local SMF derived 
by Cole et al. (2001) is also shown in Figure 4. The plot- 
ted uncertainties include Poisson errors only. The uncertain- 
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Fig. 5.— Ratio of the high-z SMFs ($) and the local SMF ($z~o.l) plot- 
ted as function of the stellar mass as measured from the maximum-likelihood 
analysis. The shaded regions represent the 1 cr uncertainties. Colors are 
as in Fig. 4. The vertical dashed and dotted lines represent the value of 
3 X 10" M0 and the z = 0.1 characteristic stellar mass, M*,^ ~ lO" Mq, 
respectively. Evidence for mass-dependent evolution is present, with the evo- 
lution to z 0.1 being larger at the low-mass end and smallest for the most 
massive galaxies. 

ties on the derived SMFs due to cosmic variance, photometric 
redshift errors, and different SED-modeling assumptions are 
quantified and discussed in § 5. 

The large surveyed area allows for the determination of the 
high-mass end with unprecedented accuracy, while the depth 
of the FIRES survey allows us to constrain also the low-mass 
end. This is particularly important because of the well-known 
correlation between the two parameters a and M*^^^. 

Figure 4 clearly shows dramatic evolution of the SMF, qual- 
itatively consistent with other studies (e.g., Fontana et al. 
2006; Perez-Gonzalez et al. 2008). The main trend is a grad- 
ual decrease with redshift of the characteristic density $*, 
rather than a change in the slope a or the characteristic stellar 
massM*^j.. The density at Mstar ^^-^ 10" has evolved by a factor 
of ~ 20 since z = 3.5, a factor of 8 since z = 2.5, and a factor 
of ^ 3.5 since z = 1.65. The data points shown in Figure 4, 
along with the best-fitting Schechter function parameters, are 
listed in Table 1 and Table 2. 

We also find evidence for mass-dependent evolution. In 
particular, the data suggest a remarkable lack of evolution 
for the most massive galaxies, with Mstar > 3 x 10", over the 
redshift range 1.3 < z < 4.0. The average density of these 
galaxies is 1.3 x 10"^ Mpc"^. The differential evolution of 
galaxies with different masses is shown more clearly in Fig- 
ure 5, where the high-redshift SMFs divided by the local SMF 
have been plotted as function of stellar mass. If the form of 
the mass function does not evolve with redshift, the curves in 
Figure 5 would be constant lines as function of stellar mass. 
On the contrary, the observed evolution of the number den- 
sities is larger for less massive galaxies and smallest for the 
most massive galaxies. 

5. UNCERTAINTIES 

The results found in the previous sections are very intrigu- 
ing. However, only Poisson errors have been considered, and, 
as previously noted, uncertainties due to photometric redshift 
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TABLE 1 

SMFS DERIVED WITH THE 1/Vmax METHOD 
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log* 





(Mpc ^ dex" 
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cr = (cTp^j + (T^j, + is the total 1 ct random error, including the Poisson errors 

(crpoi), the errors due to photometric redshift random uncertainties (crzsun', see § 5.1), 
and the error due to cosmic variance (ccv; see § 5.2); sys is the systematic uncertainty 
due to different template sets in the photometric redshift estimate, while iTsed.sys is the 
systematic uncertainty due to different SED-modeling assumptions (see § 5.1 and 5.3). 



TABLE 2 

BEST-FIT SCHECHTER FUNCTION PARAMETERS FOR THE SMFS 



Redshift 


a 






0;z,sys 


0^86(1,8X8 


(l0gA^Sar)z,sy8 


Oog^?tar)sed,sys 


^z,sys 


sea,sys 


Range 




(Mq) 


(10"^ Mpc-3 dex-1) 












z~0.1 


-1.18 ±0.03 


10.96 ±0.01 


30.87 ±4.80 














1.3 < z < 2.0 
2.0 < z < 3.0 
3.0 < z < 4.0 


„ qq+fl.13,0.23 
^■^'-0.14,0.23 
, ^,+0.30,0.50 
^■"^-0.27,0.46 
, -,Q+O.63,1.06 
^•■^^-0.55,0.98 


,„„. +0.07 ,0.11 
'"•'^-0.05.0.09 

in Qfi+o. 12:0.23 

00+0.46,1.36 


in 17+1 -7 1-3-02 
^"•^ -1.99.3.19 
, Q,+l, 44.2.43 
■^•"-1.34,2.15 
(, r-!+0.81,1.32 
"•-'-'-0.45,0.52 


+0.08 
-0.10 
+0.03 
-0.12 
+0.30 
-0.20 


+0.16 
-0.31 
+0.16 
-0.35 
+0.42 
-0.53 


+0.05 
-0.06 
+0.03 
-0.01 
+0.08 
-0.25 


+0.41 
-0.17 
+0.22 
-0.33 
+0.25 
-0.30 


+0.60 
-1.71 
+0.85 
-1.03 
+0.37 
-0.16 


+3.61 
-8.16 
+1.07 
-3.20 
+0.12 
-0.42 



The quoted errors correspond to the 1 and 2 cr errors estimated from the maximum-likelihood analysis as described in § 4.3. Also listed are 
the systematic uncertainties on the Schechter fimction parameters due to different SED-modehng assimiptions and different template sets in the 
photometric redshift estimate (see § 5.1 and 5.3). The local (z ~ 0.1) values are taken from Cole et al. (2001). 



errors (both random and systematic), cosmic variance, and 
different SED-modeling assumptions also affect the measure- 
ment of the high-redshift SMF.^-^ In this section we quantify 
the uncertainties on the measured SMFs due to these sources 
of errors, providing the first comprehensive analysis of ran- 
dom and systematic uncertainties affecting the high-z SMFs. 

5.1. Uncertainties due to photometric redshift errors 

Studies of high-redshift galaxies largely rely on photomet- 
ric redshift estimates. It is therefore important to understand 
how the photometric redshift uncertainties affect the derived 
SMFs and densities. 

We note that the results from the maximum-likelihood analysis are un- 
biased with respect to density inhomogeneities, hence not affected by cosmic 
variance. 



5.1.1. Photometric redshift random errors 

To quantify the uncertainties on the SMFs due to photo- 
metric redshift random errors we have proceeded as follows. 
First, for each galaxy in the /T-selected sample, a set of 100 
mock SEDs was created by perturbing each flux point ac- 
cording to its formal error bar. Second, we estimated pho- 
tometric redshift Zphot in the same way as described in § 2.2. 
Third, we fitted the mock SEDs to estimate stellar masses 
as described in § 3, using the default set of SED-modeUng 
assumptions. Finally, we have derived completeness limits 
in stellar mass and SMFs of galaxies with the 1 /Vmax and 
the maximum-likelihood analysis for each of the 100 Monte 
Carlo realizations of the composite ^-selected sample. This 
approach naturally addresses the fact that fainter sources tend 
to be characterized by less accurate Zphot estimates due to the 
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larger errors in their photometry, as well as sources charac- 
terized by power-law SEDs and consequently by very poorly- 
constrained Zphot estimates and very broad Zphot distributions 
derived from the Monte Carlo realizations. Moreover, the 
adopted Monte Carlo approach to estimate the uncertainties 
on the SMF due to the photometric redshift random error is to 
be preferred to the approach using the comparison of Zphot with 
Zspec, as this comparison is strongly affected by the very biased 
and incomplete sub-sample of galaxies at z > 1.5 with avail- 
able spectroscopic redshifts (see, e.g., Brammer et al. 2008). 

The contribution to the total error budget of the SMFs de- 
rived using the 1 /Vmax method due to photometric redshift 
random errors (cTz.ran) was estimated by taking, for each stel- 
lar mass bin, the lower and upper errors on <1>(M) comprising 
the central 68% of the Monte Carlo distribution. The values 
of ^ z.ran 

for each stellar mass bin in the three targeted red- 
shift intervals are listed in Table 1 . The contribution of pho- 
tometric redshift random uncertainties to the total error bud- 
get of the SMFs derived with the 1 /Vmax method is generally 
smaller (although non-negligible) than crpoi and CTcv (the er- 
ror due to cosmic variance; see § 5.2), the latter dominating 
the random error budget. This is true at all redshifts. The 
contribution of (Tz,ran is largest (although still relatively small) 
for the largest stellar mass bin (M ~ 1 1 .65), which is usually 
populated only by a handful of sources, and for the SMF at 
redshift3.0<z<4.0. 

The uncertainties on the SMFs derived using the maximum- 
likelihood analysis due to photometric redshift random errors 
appears to be negligible. This is due to the fact that, when the 
maximum-hkehhood analysis is used to derive the SMFs, the 
whole stellar mass range contributes to the determination of 
the Schechter function parameters, significantly reducing the 
impact of photometric redshift random errors. The derived 
1 a contour level from the maximum-likelihood analysis con- 
tains ~95% of the Monte Carlo realizations.'^ Therefore, the 
errors on the Schechter function parameters due to photomet- 
ric redshift random errors can be neglected. This is true for 
all three targeted redshift intervals. 

5.1.2. Photometric redshift systematic errors 

In addition to random errors, systematic errors can be 
caused by the specific choice of the templates or the tem- 
plate error function used in the estimate of the photo- 
metric redshifts. To quantify these systematic errors, we 
have repeated the Zphot estimates using the following dif- 
ferent combination of template set and template error func- 
tion: 1) eazy_vl ■O_nodust and TE.eazy_v 1 ■O_nodust, with 
eazy_vl.O_nodust equal to the default set eazy_vLO without 
the dusty template, and TE.eazy_vl .O_nodust\h& template er- 
ror function specifically constructed for the eazy_vl.O_nodust 
template set; 2) eazy_vl.O and TE.eazy_vl.O_nodust; 3) 
br07_default and TE.eazy_vl.O, with br07_default the default 
template set of Blanton & Roweis (2007). These three combi- 
nations were chosen because they resulted in Zphot-Zspec com- 
parisons of similar quahty (or only slightly worse) as that de- 
rived in § 2.2 using the default EAZY template set and tem- 
plate error function. We decided not to use the cww+kin^'^ and 
the pegaselS^^ template sets (also distributed with the EAZY 

" These results are fully consistent with the Monte Carlo simulations per- 
formed in Marchesini et al. (2007) to address the systematic uncertainties 
on the Schechter function parameters of the high-redshift rest-frame optical 
luminosity functions due to photometric redshift random uncertainties. 

The cww+kin template set comprises the empirical templates from Cole- 
man, Wu, & Weedman (1980) plus the "SBl" starburst spectrum from Kin- 



code), due to the significantly worse resulting Zphot - Zspec com- 
parisons. Modeling of the observed SEDs was then performed 
using the new sets of Zphot to derived stellar masses, and the 
completeness limits in stellar mass were then re-estimated and 
the SMFs re-derived with both the 1 /Vmax and the maximum- 
likelihood methods. The last three columns of Table 4 list 
the SMFs of galaxies at 1.3 < z < 2.0, 2.0 < z < 3.0, and 
3.0 < z < 4.0 derived using the 1 /Vmax method for each com- 
bination of template set and template error function. Table 5 
lists the best-fit Schechter function parameters a, M*^, and 
of the SMFs of galaxies at 1.3 < z < 2.0, 2.0 < z < 3.0, and 
3.0 < z < 4.0 derived using the maximum-likelihood analysis 
for each combination of template set and template error func- 
tion. 

Systematic errors were then quantified by comparing the 
resulting SMFs with the SMFs derived using the preferred 
default EAZY template set and template error function, fol- 
lowing the same approach (described in §5.3) used to quan- 
tify the systematic uncertainties due to different SED mod- 
ehng assumptions. These systematic uncertainties, C7z,sys 
for the SMFs derived with the l/Vmax method, and az,sys, 
(logM*(.„.)z.sys, and ^*^y^ for the Schechter function parame- 
ters derived with the maximum-likelihood analysis, are listed 
in Table 1 and 2, respectively. The values of CTz,sys are gen- 
erally asymmetric and larger than C7z,ran. and larger than the 
total 1 a random error a for ~l/3 of the considered stellar 
mass bins. On the contrary, the systematic uncertainties on 
the Schechter function parameters due to different template 
sets or template error functions are always smaller than the 
1 a error estimated from the maximum-likelihood analysis. 
We note that our results may be influenced by unknown sys- 
tematic effects in the redshifts, particularly at the high-mass 
end. Spectroscopic redshifts, or photometric redshifts with 
very small errors and systematics (e.g., van Dokkum et al. 
2009), are needed to confirm the shape of the mass function 
in the highest mass bins. 



5.2. Uncertainties due to cosmic variance 

As already pointed out, cosmic variance represents a sig- 
nificant source of uncertainty in deep surveys, since they are 
characterized by small areas and hence small probed volumes. 
Our composite sample is made of several independent fields 
with a large total effective area of ^ 511 arcmin^, which sig- 
nificantly reduces the uncertainties due to cosmic variance. 
Also, the large number of fields considered in this work with 
their large individual areas allows us to empirically quantify 
the field-to-field variations from one field to the other in the 
estimate of the SMF with the l/Vmax method, especially at 
the high-mass end, and to properly account for it in the error 
budget. 

In order to quantify the uncertainties due to field-to-field 
variations in the determination of the SMF, we proceeded as 
in Marchesini et al. (2007). Briefly, using the 1 /Vmax method, 
we measured <I>j, where $j is the galaxy number density in the 
stellar mass bin AM for the Jth field. For each stellar mass bin 
with n > 3, we estimated the contribution to the error budget 

ney et al. (1996); these templates have been extended in the UV and IR by 
Amouts et al. (1999) (http://www.oamp.fr/people/arnouts/LE_PHARE.html) 
The pegaselS template set contains the set of constant star formation 
rate models with additional dust reddening applied using the extinction curve 
of Calzetti et al. (2000) (see Brammer et al. 2008 for details). 
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of $ from cosmic variance using: 

rms($J) 

'7cv = (2) 

with n the number of individual fields used. For the stellar 
mass bins with « ^ 2, we adopted the mean of the rms(<l>j) 
with n ^ 3. The final 1 a random error associated to ^(M) 
is then a = (o"poi + CT^v + '''z.ran)'''^' with crpoi the Poisson error 
in each magnitude bin, and CTz ran the error due to photometric 
redshift random uncertainties as derived in § 5.1.'^ 

The values of CTcv for each stellar mass bin in the three 
targeted redshift intervals are listed in Table 1. In the red- 
shift range 1.3 < z < 2.0, cosmic variance is the dominant 
source of random errors over almost the entire probed stel- 
lar mass range, with the exception of the most massive bin 
which is populated only by a handful of sources. At z > 2.0, 
cosmic variance is generally comparable, or slightly smaller 
than, Poisson errors, due to the larger probed volumes and 
to the smaller number of galaxies with respect to the 1.3 < 
z < 2.0 redshift interval. We stress that the results from the 
maximum-likelihood analysis are not affected by cosmic vari- 
ance, since the adopted STY method is unbiased with respect 
to density inhomogeneities (e.g. Efstathiou, Ellis & Peterson 
1988). 

5.3. Systematic effects due to dfferent SED modeling 
assumptions 

As described in § 3, the default set of SED-modeling as- 
sumptions is represented by (BC03,ZQ,Kroupa,Calzetti), i.e. 
BC03 stellar population synthesis models with a pseudo- 
Kroupa (2001) IMF and solar metaUicity have been used 
in combination with the Calzetti et al. (2000) extinction 
law to derived stellar masses. With broad-band photome- 
try alone, it is not possible to constrain the metallicity, the 
IMF, the extinction law, or the stellar population synthesis 
model. Even with high-quality optical-to-MIR photometry 
and NIR spectroscopy, it is not possible to statistically con- 
strain any of the above, as shown by Muzzin et al. (2009) 
for a sample of z ~ 2 galaxies. Therefore, we have cho- 
sen (BC03,ZQ,Kroupa,Calzetti) as our default set of SED- 
modeling assumptions, instead of having the metallicity, the 
IMF, the extinction curve, and the stellar models as free pa- 
rameters. 

In the following we describe the adopted approach to quan- 
tify the systematic effects on the derived SMFs of galaxies 
due to the different choices of SED-modeling settings. 

5.3.1. Variations on the default SED-modeling assumptions 

We have derived stellar masses by fitting the observed 
SEDs with different sets of SED-modeling assumptions, by 
changing the stellar population synthesis models, the IMF, the 
metallicity, and the extinction law. 

For the additional metallicities, we have used super-solar 
(Z = 2.5 X Zq) and sub-solar (Z = 0.2 x Zq) metalUcities. 

To explore the systematic effects due to different attenua- 
tion laws, the Milky Way (MW) extinction curve by AUen 
(1976) and tiie Small Magellanic Cloud (SMC) extinction 
curve (Prevot et al. 1984; Bouchet et al. 1985) were also used. 

We note that uncertainties related to field-to-field variations could be po- 
tentially correlated with photometric redshift uncertainties. However, spec- 
troscopic complete high-mass samples at high redshift {z> 1) do not yet 
exist, resulting in very little knowledge about an existing covariance between 
density estimates and photometric redshift systematics. 



The main differences between the Calzetti et al. (2000) and 
the MW extinction laws lie in the ratio of total-to-selective ?ib- 
sorptionT^v = Av/£(£-y) (4.05 versus 3.1, respectively) and 
in the Calzetti et al. (2000) law lacking the 2175 A bump char- 
acteristic of MW dust mixtures. Otherwise, their wavelength 
dependence are fairly similar The SMC law with Ry = 2.72 
also lacks the 2175 A bump. In addition, it rises more steeply 
with decreasing wavelengths in the near-UV than the other 
two laws; in other words, the Calzetti et al. (2000) and the 
MW laws are much "grayer" at near-UV wavelengths. For 
(self-)consistency, the MW extinction law was used in com- 
bination with solar and super-solar metallicities, whereas the 
SMC curve was used with the sub-solar metalUcity. 

In addition to the pseudo-Kroupa (2001) IMF (discussed 
in § 3), we have used three additional IMFs, namely the 
Chabrier (2003) and two bottom-light IMFs.'^ Theoretical 
arguments and indirect observational evidence suggest that 
the stellar IMF may evolve with cosmic time, such that it is 
more weighted toward high-mass stars at higher redshift (see, 
e.g., Dave 2008; van Dokkum 2008; Wilkins et al. 2008). Re- 
cently, van Dokkum (2008) provided new constraints on the 
IMF at high redshift by comparing the evolution of the M /L 
ratios of early-type galaxies to their color evolution, finding a 
logarithmic slope of the IMF around 1 Mq (x = -0.3) signifi- 
cantly flatter than the present-day value (x ~ 1.3). Moreover, 
assuming a Chabrier (2003)-like parameterization of the IMF 
with an evolving characteristic mass m^, the analysis in van 
Dokkum (2008) impUes a characteristic mass Wc = 1.9 Mq at 
z = 3 - 6 (for solar metallicity). This IMF is best described as 
"bottom-light" rather than top-heavy, since it does not have 
a larger number of massive stars than a standard Chabrier 
(2003) IMF, but has a deficit of low-mass stars. For the 
bottom-light IMF, we adopted the parameterization defined 
in eq. 18 of van Dokkum (2008), with = 1.9 Mq. We also 
used a second bottom-light IMF by adopting a smaller value 
for the characteristic mass, mc = 0.3 M©. This is the charac- 
teristic mass required to reproduced the top-heavy IMF with a 
simple cutoff at 1 M© invoked by Blain et al. (1999a) for sub- 
millimeter galaxies. The Chabrier (2003) IMF is recovered by 
using OTc = 0.079 Mq. Note that the bottom-light IMFs have 
been used in combination with the Maraston (2005) stellar 
population synthesis models. 

Different stellar population synthesis models do not paint a 
consistent picture of evolution in the rest-frame NIR (probed 
by the IRAC bands). Therefore, we have explored systematic 
effects due to different stellar population synthesis models by 
performing SED modeling with the Maraston (2005) (MA05) 
and tiie Chariot & Bruzual (2008) (CB08) stellar population 
models. The BC03 and the MA05 models differ in several 
aspects: the stellar evolutionary ttacks adopted to construct 
the isochrones, the synthesis technique and the treatment of 
the thermally pulsating asymptotic giant branch (TP-AGB) 
phase. The Padova stellar tracks used in BC03 include a 
certain amount of convective-core overshooting, whereas the 
Frascati ttacks (Cassisi et al. 1997) used in MA05 do not. The 
two stellar evolutionary models also differ for the temperature 
distribution of the red giant branch phase. The differences in 
the rest-frame NIR originates mainly from a different imple- 
mentation of the TP-AGB phase (Maraston et al. 2006). Fol- 
lowing the fuel consumption approach, Maraston (2005) finds 

We note that (because of our fitting procedure) the shape of the SMF is 
identical for a Salpeter (1955) IMF, with a simple systematic shift of a factor 
of ~ 1.6 to larger stellar masses. 



The SMFs at 1 .3 < z < 4.0 and their uncertainties 
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TABLE 3 

Considered sets of SED-modeling 
assumptions 



( model,Z,IMF,dust) 

( BC03, Zq , Kroupa, Calzetti) 
(BC03,2.5 ZQ.Kroupa.Calzetti) 
(BC03,Q.2 ZQ.Kwupa. Calzetti) 
(BC03,ZQ,Kwupa,MW) 
(BC03,2.5 ZQ.Kmupa.MW) 
(BC03,Q.2 ZQ,Knmpa,SMC) 
(BC03,Zq , Chahrier, Calzetti) 
( CB08, Zq , Kroupa, Calzetti) 
(MA05,ZQ,Kroupa,Calzetti) 
(MA05,Zq, Bottom-light = 0.3,Calzetti) 
(MA05,ZQ,Bottom-light iric = l.9,Calzetti) 



The first element of the table is the default 
set of SED-modeling assumptions. 

that this phase in stellar evolution has a substantial impact on 
the NIR luminosity for ages between 0.2 and 2 Gyr. Bruzual 
& Chariot (2003) follow the isochrone synthesis approach, 
characterizing properties of the stellar population per mass 
bin. The latter method leads to smaller luminosity contribu- 
tions by TP-AGB stars. The CB08 stellar population synthesis 
models are generated with a recent version of the Bruzual & 
Chariot (2003) stellar population synthesis code which incor- 
porates a new prescription by Marigo & Girardi (2007) for 
the TP-AGB evolution of low- and intermediate-mass stars. 
Whereas the Marigo & Girardi (2007) tracks used in CB08 
account for 9 evolutionary stages in the TP-AGB (three in the 
O-rich phase, three in the C-rich phase, and three in the su- 
perwind phase), the BC03 models include only 1 evolution- 
ary stage on each of these phases. The main effect of this 
added prescription is to improve the predicted NIR colors of 
intermediate-age stellar populations (Bruzual 2007; see also 
Chariot & Bruzual 2008). 

We note that the star formation history is also a significant 
source of uncertainty. We have treated this implicitly in our 
Monte Carlo simulations, as we chose the best-fitting star for- 
mation history (out of three models) for each realization (see 
§ 3). However, it is well known that masses can be altered 
significantly by adding "maximally old" components in the 
fits and generally by allowing more complex forms of the star 
formation history than simple exponentially declining mod- 
els (e.g., Papovich et al. 2001; Wuyts et al. 2007). Fitting 
such complex star formation history models is beyond the 
scope of the present paper, but we note that multiple compo- 
nent fits tend to increase the masses, particularly for galaxies 
whose light is dominated by star bursts (see Wuyts et al. 2007 ; 
Pozzetti et al. 2007). 

The considered sets of SED-modeUng assumptions are 
summarized in Table 3. 

5.3.2. Derivation of the SMFs 

For each new combination of SED-modehng assumptions, 
we have derived the completeness Umits in stellar mass and 
the SMFs with both the 1 /Vmax method and the maximum- 
likelihood analysis. Table 4 lists the SMFs of galaxies at 
1.3 < z < 2.0, 2.0 < z < 3.0, and 3.0 < z < 4.0 derived using 
the l/Vmax method for each combination of SED-modeling 
settings. Table 5 lists the best-fit Schechter function parame- 
ters a, M*3,., and $* of the SMFs of galaxies at 1 .3 < z < 2.0, 
2.0 < z < 3.0, and 3.0 < z < 4.0 derived using the maximum- 
UkeUhood analysis for each combination of SED-modeling 



settings. In the left panel of Figures 6, the SMF of galaxies at 
1 .3 < z < 2.0 derived with the 1 /V,nax method for the default 
set (BC03,ZQ,Kroupa, Calzetti) is compared to the SMFs de- 
rived for the other considered sets of SED-modeling assump- 
tions. Similarly, the left panel of Figure 7 shows the different 
SMFs at 1.3 < z < 2.0 corresponding to the various sets of 
SED-modeling assumptions used in the maximum-UkeUhood 
analysis. 

5.3.3. The effects of dijferent SED-modeUng assumptions 

In this section we discuss in details the effects on the 
derived SMFs when changing the SED-modeling assump- 
tions. A detailed analysis of the effects of the different SED- 
modeling assumptions on the estimated stellar masses is pre- 
sented in Muzzin et al. (2009) for a sample of 34 A'-selected 
galaxies at z ^ 2. 

Stellar population synthesis models — With respect to the de- 
fault SED-modeling assumptions, using the Maraston (2005) 
models results in derived SMFs with generally steeper low- 
mass end slopes a, slightly smaller characteristic stellar 
masses M*,^^ (by < 0.1 dex), and smaller normalizations 
(by ~40%-50%). If the Chariot & Bruzual (2008) models are 
instead used, the derived SMFs have similar a, significantly 
smaller M*(2j (by ^0.1-0.2 dex), but similar $*. However, due 
to the correlation between the Schechter function parameters 
a and M*^, the SMFs derived using the Maraston (2005) and 
the Chariot & Bruzual (2008) models are in general very sim- 
ilar, resulting in a general decrease of the number densities 
of galaxies. This decrease is larger at the high-mass end, and 
smaller at the low-mass end. 

Metallicities — Changing the metalhcity from solar to sub- 
solar results in smaller characteristic densities $* by ~20%- 
30%, but no significant effect on a and M*^. Conversely, 
using super-solar metallicity results in shallower a, smaller 
*^s*tar (by ~0. 1-0.2 dex), and larger $* (by ~30%-40%). The 
SMFs derived with sub-solar metallicities are similar to those 
derived using the default SED-modehng assumptions at the 
high-mass end, but with generally smaller number densities at 
the low-mass end. The SMFs derived with super-solar metal- 
licities are instead characterized by a smaller number densi- 
ties with respect to the SMFs derived with the default SED- 
modeling assumptions. This decrease is larger at the high- 
mass end, and much smaller at the low-mass end. 

Extinction laws — Changing the adopted extinction law from 
Calzetti et al. (2000) to the MW law results in steeper a, sim- 
ilar or slightly larger M*^,^^, and significantly smaller (by 
^20%-50%). The net results on the derived SMFs is a de- 
crease in the number densities with respect to the SMFs de- 
rived with the Calzetti et al. (2000) extinction law. This de- 
crease is small at the high-mass end and much larger at the 
low-mass end. 

Using the SMC extinction curve in combination with sub- 
solar metallicity results in shghtly shallower a, smaller M*^^ 
(by ~0.1-0.15 dex), and larger (by ~40%-60%) compared 
to the SMFs derived with Calzetti et al. (2000) extinction law 
and sub-solar metallicity. The net result is a decrease of the 
number densities at the high-mass end, and an increase of 
the number densities at the low-mass end. With respect to 
the default SED-modeling assumptions, using a SMC extinc- 
tion curve in combination with sub- solar metalhcity results 
in smaller number densities at the high-mass end and similar 
number densities at the low-mass end. The latter is due to 
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Fig. 6. — Left panel: SMF of galaxies at 1.3 < z < 2.0 derived with the 1 /Vmax method. The SMF corresponding to the default set of SED-modeling assumptions 
is plotted with black filled circles and 1 cr Poisson errors; the SMFs corresponding to different sets of SED-modeling settings and different combinations of 
template sets and template error functions are plotted with different colors (no errors plotted for clarity). Right panel: SMF of galaxies at 1.3 < z < 2.0 derived 
with the 1 /Vniax method and assuming the default set of SED-modeling settings; the black error bars now include the Poisson eiTor, the error due to field-to-field 
variations, and the error due to photometric redshift random uncertainties. The gray boxes represent the total 1 a errors, with the systematic uncertainties added 
linearly to the 1 cr random errors u = (cp^j + fr^v + cr^jan) 
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Fig. 7. — Left panel: SMF of galaxies at 1.3 < z < 2.0 derived with the maximum-likelihood analysis. The SMF and its 1 cr error corresponding to the 
default set of SED-modeling assumptions is plotted with black line and gray shaded region; the SMFs corresponding to different sets of SED-modeling settings 
and different combinations of template sets and template error functions are plotted with different colors (no errors plotted for clarity); the aiTows represent the 
characteristic stellar masses M*,^. Right panel: SMF of galaxies at 1 .3 < z < 2.0 derived with the maximum-likelihood analysis and assuming the default set of 
SED-modeling settings (black solid curve); the gray shaded region represent the total 1 cr uncertainty, included the systematic uncertainties. The arrow represents 
M*!^^; the smaller eiTor bars represent the 1 cr en'or derived from the maximum-likelihood analysis; the larger error bars represent the total 1 cr eiTor,with the 
systematic uncertainties added linearly. The insert shows the parameter space (ct—M*^^^), with the best-fit values corresponding to the default set of SED-modeling 
settings (black filled circle) and its corresponding 1 and 2 cr contour levels (solid gray ellipsoids), and the best-fit values corresponding to the other SED-modeling 
assumption sets and different combinations of template sets and template error functions (colored filled circles). If the bottom-light IMFs are not considered, the 
largest systematic effects on the derived SMFs are caused by the changes in the stellar population synthesis models and the combination of super-solar metallicity 
with the MW extinction law. Much larger systematic effects are found when the bottom-light IMFs are adopted (brown and light-blue symbols), both at the high- 
and low-mass ends. 
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TABLE 4 

SMFs FROM THE 1/Vmax METHOD FOR THE DIFFERENT SED-MODELING ASSUMPTIONS 
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-5.039 
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<-5.3 
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2.0 < z < 3.0; 
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-3.309 


-3.424 


-3.451 


-3.331 


-3.354 


-3.357 


-3.417 


-3.462 


-3.741 


-3.787 


-3.368 


-3.295 


-3.309 


10.47 


-3.119 


-3.292 


-3.360 


-3.353 


-3.346 


-3.148 


-3.168 


-3.259 


-3.334 


-3.479 


-3.600 


-3.218 


-3.080 


-3.203 


10.18 


-3.401 


-3.208 


-3.188 


-3.175 


-3.210 


-3.235 


-3.229 


-3.281 


-3.131 


-3.258 


-3.534 


-3.265 


-3.435 


-3.526 


9.89 


-2.640 


-2.911 


-3.078 


-2.813 


-2.934 


-2.938 


-2.794 


-2.963 


-3.024 


-3.311 


-3.089 


-2.703 


-2.698 


-2.499 


3.0 < z < 4.0: 
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-4.809 


-4.218 


-4.273 


-4.710 


11.08 


-4.025 


-3.964 


-4.131 


-4.109 


-4.078 


-4.068 


-3.996 


-4.160 


-4.068 


-4.509 


-4.508 


-4.043 


-3.879 


-4.132 


10.79 


-3.929 


-3.886 


-4.099 


-3.778 


-3.973 


-3.814 


-3.964 


-4.012 


-3.983 


-4.073 


-4.210 


-3.903 


-3.920 


-3.892 


10.50 


-3.433 


-3.892 


-3.788 


-3.892 


-3.892 


-3.670 


-4.068 


-4.069 


-4.068 


-4.169 


-3.968 


-3.486 


-3.486 


-3.591 


10.21 


-3.141 


-3.280 


-3.092 


-2.997 


-3.070 


-3.317 


-3.214 


-3.280 


-3.070 


-4.019 


-3.547 


-3.141 


-3.141 


-3.016 



Set 1 = (BC03,ZQ,Kroupa,Calzetti), default set of SED-modeling assumptions; Set 2 = (BC03,2.5 Zq, Kwupa,Calzetti); Set 3 = 
(BC03,0.2 ZQ.Kwupa.Calzetti); Set 4 = (BC03,ZQ,Knmpa,MW); Set 5 = {BC03,2.5 ZQ.Kwupa.MW); Set 6 = (BC03,0.2 ZQ,Kroupa,SMC); 
Set 7 = (BC03,ZQ,Chabrier,Caketti); Set 8 = (CBOS.ZQ.Kroupa.Calzetti); Set 9 = (MA05,ZQ,Knmpa,Caketti); Set 10 = (MA05,ZQ,Bottom- 
light nic = 0.3,Calz.etti); Set 11 = (MA05,Zq, Bottom-light ntc = l.9,Calzetti); Set 12 = (BC03,ZQ,ZQ,Kwupa,Calzetti), with eazy_vl.O_nodust and 
TE.eazy_vl .O_nodust; Set 13 = (BC03,ZQ,ZQ,Kwupa,Calzetti), with eazy_vl.O and TE.eazy_vl.O_nodust; Set 14 = (BC03,ZQ,ZQ,Kroupa,Calzetti), 
with br07_default and TE.eazy_vl.O. 



TABLE 5 

Best- FIT Schechter parameters for the different SED-modeling assumptions 



Parametci 




Set 1 


Set 2 


Set 3 


Set 4 


Set 5 


Set 6 


Set 7 


Set 8 


Set 9 


Set 10 


Set 1 1 


Set 12 


Set 13 


Set 14 


1.3 < z < 2.0; 




































-0.99 


-0.83 


-1.05 


-1.10 


-0.96 


-0.92 


-0.94 


-1.01 


-1.17 


-1.30 


-1.24 


-1.09 


-0.99 


-0.91 


log (M*„/M0) 




10.91 


10.73 


10.97 


10.95 


10.80 


10.80 


10.84 


10.80 


10.92 


10.91 


11.32 


10.95 


10.95 


10.85 


$* [10^ 


dex-1] 


10.17 


13.78 


7.49 


7.35 


10.65 


12.72 


11.01 


9.62 


6.38 


3.64 


2.01 


8.46 


8.54 


10.77 


2.0 < z < 3.0; 




































-1.01 


-0.85 


-1.03 


-1.24 


-0.89 


-1.03 


-1.03 


-0.97 


-1.21 


-0.94 


-1.36 


-1.09 


-0.98 


-1.13 


log (M*„/M0) 




10.96 


10.83 


10.99 


10.94 


10.80 


10.88 


10.90 


10.83 


10.93 


10.62 


11.17 


10.96 


10.94 


10.98 


$* [10^ Mpc-3 


dex-1] 


3.95 


5.02 


3.14 


2.59 


4.33 


4.41 


4.10 


3.75 


2.78 


3.45 


0.75 


3.70 


4.80 


2.92 


3.0 < z < 4.0; 




































-1.39 


-1.39 


-1.31 


-1.92 


-1.74 


-1.61 


-1.49 


-0.96 


-1.44 


-1.06 


-1.69 


-1.44 


-1.09 


-1.59 


log(M*^/M0) 




11.38 


11.36 


11.36 


11.64 


11.44 


11.43 


11.41 


11.13 


11.26 


11.09 


11.41 


11.46 


11.24 


11.13 


[10"* Mpc-3 


dex-1] 


0.53 


0.42 


0.44 


0.11 


0.22 


0.34 


0.40 


0.65 


0.49 


0.42 


0.13 


0.37 


0.90 


0.65 



SED-modeling assumption sets as in Tab. 4. 
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the fact that the effects at the low-mass end of changing the 
extinction curve and the metallicity are broadly similar, but 
opposite in sign. 

IMFs~ Using the Chabrier (2003) IMF in place of the 
pseudo-Kroupa (2001) IMF does not have a significant effect 
on the derived shape of the SMFs, with only a small decrease 
in the characteristic stellar mass M*^^^. by ^^0.05 dex. 

A more complex behavior is however found when the two 
bottom-light IMFs are considered. As shown in the left 
panel of Figures 7, the shapes of the SMFs derived using 
the bottom-light IMFs are significantly different from that 
of the SMF derived with the default SED-modeling assump- 
tions. This is especially true for the bottom-light IMF with 
mc = 1.9 Mq, characterized by a steeper low-mass end and a 
characteristic stellar mass larger by a factor of 2.5. This re- 
sults is particularly important since it is commonly assumed 
that changing IMF in the SED-modeling results in a system- 
atic shift of the derived SMF, leaving the shape of the SMF 
unchanged. This is clearly not the case for bottom-light IMFs: 
the more the IMF is skewed toward high-mass stars (i.e., the 
more deficient in low-mass stars the IMF is), the larger is the 
effect on the shape of the derived SMF. 

Another very interesting result is the resulting higher num- 
ber density of massive galaxies when using the bottom-light 
IMF with mc =1.9 M© with respect to the SMFs derived with 
the other IMFs. This result might come unexpectedly at first. 
Naively, one would expect that, by making the IMF more de- 
ficient in low-mass stars, which dominate the stellar mass of 
a galaxy but contribute little to the integrated light, the de- 
rived stellar masses would be smaller compared to those de- 
rived from the other considered IMFs, as a consequence of 
a lowering of the M/L ratio. However, as already pointed 
out by van Dokkum (2008), the number of turn-off stars is 
also reduced for > 0.4 Mq, and these stars dominate the 
light at rest-frame optical wavelengths. Moreover, the turn- 
off mass can be similar to m^, meaning that the effect on 
the M/L ratio is not a constant, but depends on the age of 
the population. A final complication is the mass in stellar 
remnants, which is a larger fraction of the total stellar mass 
for more top-heavy IMFs. Using simple stellar evolutionary 
tracks but not full stellar population synthesis modeling, van 
Dokkum (2008) calculated the effects of changing character- 
istic mass on the M/Ly ratio for stellar populations of dif- 
ferent ages, from 0.1 to 10 Gyr They found that for young 
ages the M/L ratio steadily declines with increasing mc, but 
the behavior is more complex when mc becomes similar to the 
turn-off mass. Specifically, they found that for mc ^ 1 Mq and 
old ages the mass function becomes remnant dominated, and 
the M/L ratios approach, or even exceed, those implied by a 
Salpeter (1955) IMF. We can directly test their conclusions 
by correctly treating the above issues with the available stel- 
lar population synthesis models constructed with the bottom- 
light IMFs. The effect of changing characteristic mass on the 
M/Ly ratio for different population ages is shown in Figure 8 . 

These results are consistent with those obtained by van 
Dokkum (2008), confirming that the M/Ly ratios for old stel- 
lar populations and high characteristic mass can approach, 
and even exceed, the M/Ly ratios implied by Chabrier (2003) 
and Salpeter (1955) IMFs, due to the stellar population be- 
coming remnant-dominated. 

The SMFs derived with the bottom-light IMFs in the SED- 
modeling assumptions can be now easily explained with the 
shown behavior of the M/Ly ratio in mind. When an IMF 
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Fig. 8. — Effect of changing characteristic mass on the M/L\ ratio for 
stellar populations of ages 0.1, 0.3, 0.5, 3, 5, and 10 Gyr (from bottom to 
top) using stellar population synthesis models (filled symbols), following van 
Dokkum (2008). The solid curves, derived using simple stellar evolutional^ 
tracks, were taken from van Dokkum (2008). Three different chai'acteristic 
masses have been considered: nic = 0.08 (i.e., Chabrier 2003 IMF; dark green 
circle), mc = 0.3 (brown triangles), and = 1.9 Mq (blue squares). For 
mc = 1.9 Mq and old ages, the mass function becomes remnant dominated, 
and the M/L ratios approach, and even exceed, those implied by a Salpeter 
(1955) IMF (gray solid line). 

with mc = 0.3 is adopted, the M/L ratios are always smaller, 
or at most comparable, to the M/L ratios derived with a 
Chabrier-like IMFs. Therefore, the derived stellar masses are 
always smaller, and the derived SMF is, broadly speaking, 
shifted to smaller masses. Note, however, that the depen- 
dence of the M/L ratio on age will affect the specific shape 
of the SMF. When an IMF with mc = 1 .9 is adopted, the M/L 
ratios are larger than those derived with a Chabrier (2003) 
IMF when the age of the population is larger than r^O.9 Gyr. 
Consequently, those galaxies which are fitted by ages older 
than 0.9 Gyr will have much larger stellar masses than those 
derived assuming a Chabrier (2003) IMF. On the contrary, 
those galaxies fitted by ages younger than r^O.9 Gyr will have 
smaller stellar masses compared to those derived assuming a 
Chabrier (2003) IMF The net effect on the derived SMF is a 
significant increase of the number densities of massive galax- 
ies, which are usually characterized by old stellar population, 
and a decrease of the number densities of low-mass galaxies, 
which are usually characterized by young stellar population. 

Note that the net effect on the derived SMF caused by as- 
suming a bottom-light IMF is also a function of redshift, since 
the maximal age of the stellar population is limited by the age 
of the universe at that redshift. As the age of the stellar popu- 
lations get younger going to higher redshifts, the effect on the 
SMF due to a a bottom-light IMF will be closer to a system- 
atic shift to smaller stellar masses without changing the shape 
significantly. 

Summary — Broadly speaking, the different combinations of 
SED-modeling assumptions result in smaller estimates of the 
stellar masses with respect to the stellar masses derived using 
the default set. Consequently, the systematic effects on the 
SMFs are largest at the high-mass end of the SMFs, due to its 
steep slope and rapid changes in number density as function of 
stellar mass. The net effect on the derived SMFs is an average 
decrease of the number densities of galaxies at the high-mass 
end, while the systematic effects are generally smaller at the 
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low-mass end. If the bottom-light IMFs are not considered, 
the largest systematic effects are caused by the changes in 
the stellar population synthesis models and the combination 
of super-solar metallicity with the MW extinction law. The 
largest systematic effects are caused by the use of the bottom- 
light IMFs. 

5.3.4. The systematic uncertainties in the SMF due to 

different SED-modeling assumptions 

The systematic effects on the SMFs due to different SED- 
modeling assumptions have been quantified by comparing the 
resulting SMFs with those derived using the default set of set- 
tings (BC03,ZQ,Kroupa,Calzetti). Note that we implicitly as- 
sume that changes in the derived SMFs are the result of the 
changes we made to the model parameters. We cannot ex- 
clude subtle second-order effects that may influence the fit- 
ting procedure, but given the excellent agreement between 
the maximum-likelihood and the 1 /Vmax estimators these are 
likely much smaUer than the effects that we are measuring 
here. 

For the l/Vmax method, the systematic uncertainties of 
$(M) have been estimated by taking, for each stellar mass 
bin, the difference between the maximum (and minimum) 
value of <I>(M) allowed by all the considered combinations of 
SED-modehng settings and the value of $(M) derived with 
the default set. These systematic uncertainties (tTsys) are Usted 
in Table 1 and were then added linearly to the 1 a errors 
cr = (cpoi + CTcv + cr^rjin)'^^ (which include the Poisson error, the 
error due to field-to-field variations, and the error due to pho- 
tometric redshift random uncertainties) to obtain the total 1 a 
errors. In the right panel of Figure 6, we show the SMF of 
galaxies at 1.3 < z < 2.0, plotting the 1 a errors with and 
without the contribution of the systematic effects due to dif- 
ferent SED-modeling assumptions and different combinations 
of template sets and template error functions. 

For the maximum-UkeUhood analysis, the systematic un- 
certainties on the Schechter function parameters have been 
estimated by taking the difference between the maximum and 
minimum values derived when using all the considered com- 
binations of SED-modehng settings and the value correspond- 
ing to the default set. These systematic uncertainties (asys, 
M*yj, and $*ys), are listed in Table 2. The right panel of Fig- 
ure 7 shows the SMF of galaxies at 1 .3 < z < 2.0 derived with 
the default set of SED-modeling settings and the total 1 ct un- 
certainties after including the systematic uncertainties due to 
different SED-modeling assumptions and different combina- 
tions of template sets and template error functions; also plot- 
ted is the parameter space ia-M*^^). 

5.4. Stellar Mass Functions with all uncertainties 

Figure 9 shows the evolution of the SMF of galaxies from 
z = 4.0 to z = 1.3 including the contribution of random and 
systematic uncertainties in the error budget, i.e., the Poisson 
errors, the uncertainties due to cosmic variance and photo- 
metric redshift random errors, and the systematic uncertain- 
ties due to different SED-modeling assumptions and different 
combinations of template sets and template error functions. 
These errors are also listed in Table 1 and 2. Most of the sys- 
tematic effects are in the same direction, with a resulting net 
effect of decreasing the observed number densities, especially 
at the high-mass end and in the highest targeted redshift inter- 
val. The (M*^- a) plane is also plotted in Figure 10 showing 
the effect of systematic uncertainties on the Schechter func- 
tion parameters. 



The systematics uncertainties are a dominant contribution 
to the overall error budget. The largest contribution to the sys- 
tematic uncertainties due to different SED-modeling assump- 
tions are due to changes in the adopted IMF, specifically when 
a bottom-light IMF is used. The systematic uncertainties due 
to different combinations of template sets and template error 
functions in the estimate of the photometric redshifts are al- 
ways smaller than the systematic uncertainties due to different 
SED-modeling assumptions, especially when the maximum- 
UkeUhood analysis is used. The maximum-likelihood analysis 
is indeed quite robust against photometric redshift errors, both 
random and systematic, and the dominant source of uncertain- 
ties are the systematic errors due to different SED-modeling 
assumptions. This is true at all redshifts but the highest red- 
shift range, where Poisson errors represent a significant con- 
tribution to the error budget. As shown in the insert of Fig- 
ure 7 for the redshift range 1.3 < z < 2.0, the changes in 
the Schechter function parameters when using different SED- 
modeling assumptions, in comparison with the random errors, 
are very significant (> 2 a). At 2.0 < z < 3.0, the changes 
are sUghtly less significant, but still at the ~ 2 tr level, while 
at 3.0 < z < 4.0, where Poisson uncertainties are very large, 
the changes are mostly at the 1 a level. When the 1 /Vmnx 
method is used, cosmic variance is the dominant source of 
random errors at 1.3 < z < 2.0 in all stellar mass bins, but 
it becomes comparable to the Poisson errors at 2.0 < z < 3.0. 
The contribution of photometric redshift random uncertainties 
to the total error budget is generally smaller than Poisson er- 
rors, and increases going to higher redshifts. The relative con- 
tribution of systematic uncertainties is smallest at the highest 
targeted redshift interval, 3.0 < z < 4.0, where random errors 
contribute significantly to the total error budget. 

If the systematic uncertainties are included, the results high- 
lighted in § 4.4 are no longer robust. In particular, we cannot 
exclude a strong evolution (by as much as a factor of ~ 50) in 
the number density of the most massive (Mstar > 10^'-^) galax- 
ies from z = 4.0toz = 1.3. We note that the effects of system- 
atic uncertainties due to different SED-modeling assumptions 
are likely smaller when the redshift evolution is considered, 
as some errors would cancel out when comparing the SMFs 
at two different epochs. However, as the metallicity, IMF, and 
appropriate extinction law may all evolve with redshift, it is 
unclear to what extent this cancellation of errors actually oc- 
curs. 

5.5. Comparison with previous works 

Figure 1 1 shows the comparison of the SMFs derived in this 
study with other works from the literature (see Appendix C for 
a detailed discussion of the individual works). As discussed 
above, cosmic variance and systematic uncertainties dominate 
the errors. However, most literature studies do not give esti- 
mates for these errors. Therefore, we show the comparison 
twice in Figure 1 1 : the panels in the top two rows include 
only Poisson errors and uncertainties due to photometric red- 
shift random errors, while the panels in the bottom two rows 
include all sources of error (with the exclusion of the sys- 
tematic effects due to the bottom-light IMFs). To highlight 
similarities and differences between our SMFs and those de- 
rived in other works, we also plot A$ = log $others - log $ours 
as function of stellar mass in the second- and fourth-row pan- 
els. We note that the error bars of the different surveys cannot 
be compared directly, as they were not derived in a uniform 
way. 

From the top panels of Figure 11, it is obvious that our 
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Fig. 9. — SMFs of galaxies at redshift 1.3 < z < 2.0 (blue), 2.0 < z < 3.0 (green), and 3.0 < z < 4.0 (red). Left panel: SMFs of galaxies derived using the 
1 /Vmax method (filled circles); the error bars include Poisson errors, photometric redshift uncertainties, and en'ors due to cosmic variance. The shaded boxes 
(orange, green, and cyan, corresponding to the redshift intervals 3.0 < z < 4.0, 2.0 < z < 3.0, and 1.3 < z < 2.0, respectively) represent the total 1 a uncertainties 
in the measurements of the SMFs as described in § 5, with systematic errors added linearly to the plotted error bars. Right panel: SMFs of galaxies derived using 
the maximum-likelihood analysis (solid curves); the shaded regions represent the total 1 a uncertainties as described in § 5, including the systematic uncertainties. 
The arrows show the best estimates of M^j^,., with their en'or bars including also systematic uncertainties. 

from the literature, and also between the works in the liter- 
ature themselves, for some redshift and stellar mass ranges. 
The disagreements between the different SMFs increase with 
increasing redshift. Our SMFs fall somewhere in the middle 
of the SMFs from the literature. The largest disagreement is 
with the SMFs from Drory et al. (2005) at the low-mass end at 
all redshifts. At the high-mass end, the largest disagreement 
is with the SMFs from Fontana et al. (2006) at z ~ 3.5, from 
Perez-Gonzalez et al. (2008) at z ^ 2.5, and from Eisner et al. 
(2008) at z ^ 1.6. The large differences between the SMFs 
from Fontana et al. (2006) and Eisner et al. (2008) are in- 
teresting, since both were derived from the GOODS-MUSIC 
catalog. The former was derived from a /T-selected catalog, 
while the latter from a z-selected catalog. We note, however, 
that Fontana et al. (2006) claim that their z-selected SMF is 
very similar to their /^-selected one. 

Once the systematic uncertainties are taken into account, as 
shown in the bottom panels of Figure 11, the SMFs derived 
in this work become consistent with most of the SMFs from 
the literature. The low-mass end of the SMFs at z ^ 2.5 and 
z ^ 1.6 from Drory et al. (2005) is still significantly steeper 
with respect to both our SMFs and the other SMFs from the 
literature. A possible explanation is the different way the 
completeness limits in stellar mass has been derived by Drory 
et al. (2005) (SSP-derived completeness), as this potentially 
over-corrects densities at the low-mass end. 

We stress that most of the disagreements between the dif- 
ferent measurements of the SMFs stem from an incomplete 
analysis of the errors. The errors on the SMFs from the lit- 
erature include only Poisson errors (e.g., Drory et al. 2005), 
or Poisson errors and errors from photometric redshift uncer- 
tainties (but not cosmic variance, e.g. Fontana et al. 2006; 
Perez-Gonzalez et al. 2008; Eisner et al. 2008). Field-to-field 
variations is a significant source of errors when the SMF is 
derived using the 1 /Vmax method. This is true at all redshifts. 
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Fig. 10. — Parameter space (a-M*,;,,.) derived from the maximum- 
likelihood analysis. The red, dark green, and blue filled circles are the 
best-fit values of a and M*^.^ at redshift 3.0 < z < 4.0, 2.0 < z < 3.0, and 
1.3 < z < 2.0, respectively. The red, dark green, and blue curves represent 
their 1 and 2 cr contour levels, respectively. The filled regions show the 1 a 
allowed values of a and M*^.^^ after inclusion of the systematic uncertainties 
in the error analysis. The black filled square represents the redshift z ~ 0. 1 
value from Cole et al. (2001). 



SMFs agree with those from the literature for some redshift 
and stellar mass ranges, but disagree for others. Our SMFs 
are generally in good agreement with those from Eisner et 
al. (2008), Perez-Gonzalez et al. (2008), and Pozzetti et al. 
(2007). Broad agreement is also found with the SMFs at z < 3 
from Fontana et al. (2006), and with the high-mass end of the 
SMFs at z < 3 from Drory et al. (2005). However, there is 
also significant disagreement between our SMFs and those 
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11 10 11 10 11 10 

Fig. 1 1. — Comparison between the SMFs derived from this work and previous measurements from the literature. First row: The SMFs derived from this 
work are shown as filled red, dark green, and blue circles (1/Vmax method) and solid curves (maximum-likelihood analysis). The 1 a error bars of the l/Vmax 
measurements include Poisson errors and uncertainties from photometric redshift random uncertainties, but not cosmic variance and systematic uncertainties. 
Similarly, the 1 cr error of the maximum-likelihood measurements (orange, green, and cyan shaded regions) do not include systematic uncertainties. Previous 
works are plotted as filled stars and dashed curves (Fontana et al. 2006; F06); open circles and solid curves (Perez-Gonzalez et al. 2008; P08); open stars and 
dot-dashed curves (Eisner et al. 2008; E08); open triangles and dotted curves (Drory et al. 2005; D05); open squares and long dot-dashed curves (Pozzetti et 
al. 2007; P07). Second row: Symbols as in the first row panels, but now the differences between the SMFs from the literature and those derived in this work, 
A$ = log<I>o,hers-log'I> ouisi plotted as a function of stellar mass. Third and forth row panels: Symbols as in first and second row panels, respectively, with 
cosmic variance added in quadrature to the error bars, and systematic uncertainties now included in the total error budget represented by the shaded gray boxes (for 
the 1 /Vmax points) and shaded orange, green, and cyan regions (for the maximum-likelihood measurements). The systematic effects due to the bottom-light IMFs 
are not included. Most of the disagreements between the different measurements of the SMFs stem from an incomplete analysis of the errors. A comprehensive 
analysis of random and systematic uncertainties is necessary to reconcile the different measurements of the high-z SMFs. 
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especially at the high-mass end, but cosmic variance domi- 
nates the error budget at z '--^ 1.6. The maximum-likelihood 
estimator, unbiased with respect to density inhomogeneities, 
has been applied only by Fontana et al. (2006) and Pozzetti et 
al. (2007), while the other works have simply fitted the SMFs 
derived with the 1 /Vmax with a Schechter function. Finally, it 
is extremely important to include the systematic uncertainties 
due to different SED-modeling assumptions, which dominate 
the total error budget and are necessary to reconcile the dif- 
ferent measurements of the high-z SMFs. 

6. DENSITIES 

In this section we present estimates of the stellar mass den- 
sity pstar derived by integrating the best-fit Schechter function 
obtained from the maximum-likelihood analysis (no use of 
the l/Vmax results has been done in the estimate of the stel- 
lar mass densities). The stellar mass density (obtained by 
integrating the SMF derived from the maximum-likelihood 
analysis) is a more robust measurement than the Schechter 
function parameters M*;^^, a, and <!>*, because the errors in 
these parameters are highly correlated. Stellar mass densities 
have been estimated adopting three different integration inter- 
vals: Pstar^^'"', where the integration was performed for stel- 
lar masses 10^ < M^i^^/Mq < 10'^; pl^'^"^^\ where the inte- 
gration was performed for stellar masses 10'" < M^i^^/Mq < 



10''; and pi 



11<M<12 



, where the integration was performed for 



stellar masses 10" <Mstar/A^0 < 10'^, i.e. massive galaxies. 
These estimates are listed in Table 6, along with the 1 a errors 
and the values of the stellar mass density at z ~ 0.1 estimated 
from the SMF of Cole et al. (2001). Note that the contribu- 
tion of galaxies less massive thanMstar = 10*^ M0 to the global 
stellar mass density is negligible if the Schechter parameteri- 
zation of the SMF is a good approximation and valid also at 
stellar masses smaller than probed by our composite sample. 

The 1 a errors on the stellar mass densities have been esti- 
mated by deriving the distribution of all of the values of pstar 
allowed within the 1 a solutions of the Schechter function 
parameters from the maximum-likelihood analysis. The con- 
tribution to the total error budget from photometric redshift 
random uncertainties (derived from the 100 Monte Carlo real- 
izations described in § 5.1) was added in quadrature. The 1 a 
errors including the systematic uncertainties were estimated 
in the same way by deriving the distribution of all of the val- 
ues of pstai- allowed within the 1 <j solutions obtained using 
different SED-modeling assumptions and different combina- 
tions of template sets and template error functions to estimate 
the photometric redshifts. 

In Figure 12 we show the evolution of the total stellar mass 
density as a function of redshift, together with a compilation 
of results from the literature. The values from the literature, 
derived assuming a Salpeter (1955) IMF, have been scaled to a 
pseudo-Kroupa (2001) IMF by dividing the stellar mass den- 
sities by a factor of 1.6. The values from the literature of 
the stellar mass density converted to our IMF are listed in Ta- 
ble 7. The measured evolution of the total stellar mass density 
from z = 4.0 to z = 1 .3 is broadly consistent with most previ- 
ous measurements in the literature. Our measurements are 
among the currently most accurate measurements of the total 
stellar mass density at these redshifts. This is due to the large 
surveyed area, the large number of independent fields, and 
the high-quality of the optical-to-MlR data. Only the value at 
z ^ 3.5 is characterized by a large random error, due to the 
large uncertainties on the low-mass end slope of the derived 
SMF. Moreover, our measurements are the first to include a 
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Fig. 12. — Evolution as a function of redshift of the global stellar 
mass density estimated by integrating the SMFs over the stellar mass range 
10** < Mscji/Mq < 10'^ (the contribution to the total stellar mass density of 
galaxies with Mstar < 10* Mq is negligible if the Schechter parameteriza- 
tion of the SMF is valid also at stellar masses smaller than probed by our 
sample). Red symbols represent the total stellar mass densities estimated in 
this work (shaded boxes do not include the systematic uncertainties while er- 
ror bars do; the thick error bars do not include the systematic effect of the 
two bottom-light IMFs, while the thin error bars do). The estimates of the 
total stellar mass densities from the literature were taken from Cole et al. 
(2001) (open square); Dickinson et al. (2003) (filled green triangles; thick er- 
ror bars include random errors and cosmic variance alone, while the thin error 
bars include the systematic uncertainties due to different SED-modeling as- 
sumptions); Glazebrook et al. (2004) (open upside-down triangles; error bars 
including Poisson errors only); Drory et al. (2005) (open circles); Eisner et 
al. (2008) (filled cyan circles); Fontana et al. (2006) (filled circles); Perez- 
Gonzalez et al. (2008) (filled dark green squares); and Rudnick et al. (2006) 
(filled upside-down triangles; error bars including random errors and uncer- 
tainties due to cosmic variance). For the measurements from the literature, 
only Poisson errors and errors due to photometric redshift uncertainties (as 
derived in the corresponding work), are plotted without uncertainties due to 
cosmic variance and systematic errors, unless stated otherwise. The horizon- 
tal dotted lines represent 50%, 25%, 10%, 5%, and 2% (from top to bottom, 
respectively) of the total stellar mass density at z = 0.1. The stellar mass 
density has increased by a factor ~17;';jq, ~9itl, and ~4it0.2 from z = 3.5, 
z = 2.5, and z = 1.65, respectively, down to z = 0.1. A much stronger evo- 
lution with redshift of the stellar mass density is however allowed once the 
systematic uncertainties are taken into account. 

comprehensive analysis of random and systematic uncertain- 
ties in the error budget. Note that the mass density at z > 3 
is very poorly constrained when systematic uncertainties are 
included. We stress that this arises despite the fact that our 
sample is arguably the best suited for studying the total stellar 
mass density (as it samples both the low-mass and high-mass 
end in a homogeneous way). All previous studies in the liter- 
ature suffer from similar, or larger, uncertainties. 

The total stellar mass density has increased by a factor of 
- 17!^0' 9 ± 1, and ~ 4 ± 0.2, from redshift z = 3.5, z = 2.5, 
and z = 1.65, respectively, down to redshift z ^ 0.1. Due to 
the systematic uncertainties on the derived stellar mass den- 
sities (represented with thin error bars in Figure 12), a much 
stronger evolution with redshift of the total stellar mass den- 
sity than previously measured is actually allowed. Systematic 
uncertainties still allow for an increase as large as a factor of 
~ 70, ~ 25, and ^ 9 of the global stellar mass density (^ 44, 
~ 14, and ~ 6, if the effects of the bottom-light IMFs are 
excluded) from redshift z = 3.5, z = 2.5, and z = 1 .65, respec- 
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TABLE 6 
Stellar mass densities 



Redshift 


logp|<M<'3 


Iogpio<M<" 


logpii,<M<i^ 


Range 


(Mq 


(M© Mpc-3) 


(M0 Mpc-3) 


z~0.1 


8.51 ±0.07 


8.27 ±0.07 


7.92 ±0.07 



1.3 <z<2.0 
2.0<z<3.0 
3.0<z<4.0 



7.91 
7.55 
7.27 



+0.02(0.02.0.02) 
0.02(0.35.0.16) 
+0.0.5(0,12.0.12) 
-0.04(0.43.0.18) 
.+0.37(0.93,0.93) 
-0.13(0.61,0.39) 



7. 

7.31 
6, 



+0.02(0.04.0.04) 
-0.02(0.40.0.16) 
+0.05(0.11.0.11) 
-0.05(0.48.0.14) 
,+0.19(0.19;0.19) 
'-0.22(0.67,0.50) 



7.38: 
7.07 
6.91 



+0.07(0.09.0.08) 
-0.06(0.62,0.38) 
+0.08(0.13.0.13) 
-0.09(1.10.0.39) 
+0.13(0.13,0.13) 
-0.16(0.94,0.73) 



Stellar mass density estimated by integrating the best-fit Schechter 
SMF over the specified stellar mass range. The quoted 1 a errors 
include Poisson errors and errors due to photometric redshift uncer- 
tainties; the numbers in parenthesis are the 1 a errors including the 
systematic uncertainties due to different SED-modeling assumptions 
and different combination of template set and template error function 
used in the estimate of the photometric redshifts; in the second num- 
ber in parenthesis, the effect of the bottom-light IMFs is excluded. 



TABLE 7 

Stellar mass densities from the literature 



Redshift 




Redshift 




Range 


(Mq Mpc-3) 


Range 


(Mq Mpc-3) 





Dicldnson et al. (2003) 


Glazebroolc et al. (2004) 


0.6 


<z< 1.4 


8.26±0.08 


0.8 


<z< 1.1 


'•^°-0.14 


1.4 


< z < 2.0 


^ g^+O.17(0.33) 
'■°"-0.13(0.22) 


1.1 


<z< 1.3 


'■"^-0.14 


2.0 


<z< 2.5 


7 i-o+0. 11(0.54) 
'■-^ -0.07(0.29) 


1.3 


< z < 1.6 


7 gn+0.14 
'■^"-0.14 


2.5 


< z < 3.0 


7 -.,+0.23(0.48) 
-0.14(0.51) 


1.6 


< z < 2.0 


7 4Q+0.12 
'■^^-0.14 




Drory et al. 


(2005) 




Fontana et al. (2006) 


0.25 


< z < 0.75 


8.30±0.15 


0.4 


< z < 0.6 


8.26±0.03 


0.75 


<z< 1.25 


8.16±0.15 


0.6 


< z < 0.8 


8.17±0.02 


1.25 


<z< 1.75 


8.00±0.16 


0.8 


<z< 1.0 


8.09±0.03 


1.75 


<z< 2.25 


7.85±0.20 


1.0 


<z< 1.3 


7.98±0.02 


2.25 


< z < 3.00 


7.75±0.20 


1.3 


<z< 1.6 


7.87±0.05 


3.00 


< z < 4.00 


7.58±0.20 


1.6 


< z < 2.0 


7.74±0.04 


P^rez-GonzSlez 


et al. (2008) 


2.0 


<z<3.0 


7.48±0.04 


0.2 


< z < 0.4 


8.41 ±0.06 


3.0 


< z < 4.0 


'■"'-0.15 


0.4 


<z<0.6 


8.37±0.04 




Rudnick et al. (2006) 


0.6 


<z<0.8 


8.32±0.05 


0.0 


<z< 1.0 


o 17+0.05 

7 99+0.05 


0.8 


<z< 1.0 


8.24±0.05 


1.0 


< z < 1.6 


1.0 


<z< 1.3 


8.15±0.05 


1.6 


< z < 2.4 


7 oo+8:3^ 


1.3 


< z < 1.6 


7.95±0.07 


2.4 


<z< 3.2 


7 7i+0.()8 
'■'^-0.43 


1.6 


< z < 2.0 


7.82±0.07 




Eisner et al. 


(2008) 


2.0 


<z<2.5 


7.67±0.08 


0.25 


< z < 0.75 


8.37±0.03 


2.5 


<z<3.0 


7.56±0.18 


0.75 


<z< 1.25 


8.17±0.02 


3.0 


<z<3.5 


7.43±0.14 


1.25 


<z< 1.75 


8.02±0.03 


3.5 


< z < 4.0 


7.29±0.13 


1.75 
2.25 
3.00 


< z < 2.25 

< z < 3.00 

< z < 4.00 


7.90±0.04 
7.73±0.04 
7.39±0.05 



All the values from the literature have been divided by 1.6 to 
convert to the pseudo-Kroupa (2001) IMF used in this work. 
The quoted errors include only Poisson errors and errors due to 
photometric redsliift uncertainties, without uncertainties due to 
cosmic variance and systematic errors, unless stated otherwise. 
The quoted errors of Dickinson et al. (2003) include random 
errors and cosmic variance, with the numbers in parenthesis 
including the systematic uncertainties due to different SED- 
modeling assumptions; The quoted errors of Rudnick et al. 
(2006) include random errors and uncertainties due to cosmic 
variance. 

lively, to z ^ 0.1. However, the effects of systematic uncer- 
tainties due to different SED-modeling assumptions are likely 
smaller when the redshift evolution is considered, as some 
errors would cancel out when comparing the stellar mass den- 
sities at two different epochs. 



If only random errors are taken into account, there are a 
few measurements from the literature that are in significant 
disagreement with ours. The stellar mass density in Eisner et 
al. (2008) at z ^ 2.5 is significantly larger than our measure- 
ment by about 60%. We note however that the error bars of the 
measurements from Eisner et al. (2008) only include random 
errors, but not uncertainties from cosmic variance and sys- 
tematic errors. Because the analysis of Eisner et al. (2008) is 
based on the single, relatively small, field of GOODS-CDFS, 
field-to-field variance is a particular important source of er- 
rors. We stress the importance of a comprehensive analysis 
of the errors, including systematic uncertainties. Excluded 
our works, the only measurements of the stellar mass density 
including a complete analysis of all the errors, including sys- 
tematic uncertainties, come from the work of Dickinson et al. 
(2003). Because Dickinson et al. (2003) include both cosmic 
variance and systematic uncertainties in their error budget, 
they are consistent with ours, despite their estimates of the 
stellar mass densities lying systematically below ours. Our 
measurements are characterized by much smaller errors due 
to our much larger surveyed area (by a factor of ~ 80). 

We finally note that at z > 2, we cannot constrain the low- 
mass end of the SMFs very well, becoming incomplete at 
~ lO'" Mq. If the SMF is much steeper at the low-mass 
end than indicated by our maximum-likelihood analysis, we 
could be missing a significant fraction of the integrated stel- 
lar mass. Very recently, Reddy & Steidel (2009) suggested 
that up to ~50% of the total stellar mass at 1 .9 < z < 3.4 is in 
faint-UV galaxies with masses smaller than ^ 10'" Mr;) (com- 
pared to ^10%-20% from an extrapolation of our Schechter 
fits). However, Reddy & Steidel (2009) do not measure stel- 
lar masses directly, but convert UV luminosity to stellar mass. 
Significantly deeper NIR observations are required to directly 
probe and constrain the low-mass end of the high-z SMF. 

7. COMPARISON WTTH MODEL PREDICTIONS 

In this section we compared the derived SMFs of galaxies 
at 1 .3 < z < 4.0 with those predicted by the latest generation 
of galaxy formation models. Specifically, we have included 
the predictions from the semi-analytic models of Monaco et 
al. (2007), SomerviUe et al. (2008), and Wang et al. (2008), 
which include active galactic nuclei (AGN) feedback. We re- 
fer to those papers for detailed descriptions of their models, 
and to Fontanot et al. (2009) for a detailed comparison of 
these models. A Chabrier (2003) IMF has been assumed in all 
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models. We have therefore scaled their predictions to match 
our pseudo-Kroupa IMF by multiplying their stellar masses 
by 1.12. The model-predicted SMFs have been convolved 
with a normal distribution of standard deviation 0.25 dex, in- 
tended to represent measurement errors in logMstar- 

The semi-analytic model of Monaco et al. (2007), MOR- 
GANA, attempts, through modeling of cooling, star forma- 
tion, feedback, galactic winds and superwinds, AGN activity 
and AGN feedback, to move from a phenomenological de- 
scription of galaxy formation to a fully physically motivated 
one. We refer to Monaco et al. (2007) for a detailed descrip- 
tion of all physical processes included in MORGANA. The 
predicted SMFs and stellar mass densities from the MOR- 
GANA model adopted in this work were derived assuming 
a WMAP-3 cosmology (Spergel et al. 2007) and including 
some minor improvements with respect to Monaco et al. 
(2007) (Lo Faro et al., in prep.). 

The model predictions from Wang et al. (2008) were de- 
rived using the Garching semi-analytic model implemented 
on the Millennium dark matter simulation described in 
Springel et al. (2005). Specifically, the semi-analytic model 
described in De Lucia & Blaizot (2007) was used, which built 
on previous works by the "Munich" galaxy formation group 
(see Kauffmann & Heahnelt 2000; Springel et al. 2001; and 
De Lucia et al. 2004 for detailed descriptions of the scheme 
for building the merger tree and the prescriptions adopted 
to model the baryonic physics, most notably those associ- 
ated with the growth of and the feedback from black holes in 
galaxy nuclei and the cooling model). Here, we specifically 
consider the "C" model in Wang et al. (2008), which assumes 
a WMAP-3 cosmology. This change in cosmology results in 
a significant delay of structure formation in comparison with 
WMAP-1 results. Therefore, to compensate for the delay in 
structure formation, model "C" has twice as much the star for- 
mation efficiency with respect to the efficiency assumed by 
De Lucia & Blaizot (2007). This increase in efficiency has to 
be compensated by much higher feedback efficiencies (both 
from supernovae and from AGN) to prevent the overproduc- 
tion of stars at late times. 

The semi-analytic model of Somerville et al. (2008), built 
on the previous models described in Somerville & Primack 
(1999) and Somerville et al. (2001), presents several improve- 
ments, including, but not limited to, tracking of a diffuse stel- 
lar halo component built up of tidally destroyed satellites and 
stars scattered in mergers, galaxy-scale AGN-driven winds, 
fueling of black holes with hot gas via Bondi accretion, and 
heating by radio jets. The prediction from the Somerville et 
al. (2008) semi-analytic model are taken from their fiducial 
model WMAP-3 model, which adopts a fraction /scatter = 0.4 
of the stars in merged satellite galaxies added to a diffuse 
component distributed in a very extended halo or envelope. 

Figures 13 and 14 show the comparison between the SMFs 
measured in this work (with the exclusion of the systematic 
effects due to the bottom-light IMFs) and those predicted by 
the models. Figure 13 shows the comparison for the model- 
predicted SMFs without convolution with a normal distribu- 
tion, while Figure 14 shows the comparison for the model- 
predicted SMFs convolved with a normal distribution of stan- 
dard deviation 0.25 dex. To highlight similarities and dif- 
ferences between our SMFs and the model-predicted SMFs, 
we also plot A$ = log$modeis-log$ours as function of stel- 
lar mass in the bottom panels of Figures 13 and 14. Broadly 
speaking, the predicted SMFs are too steep with respect to 
the observed SMFs. If the comparison is done with the not- 



convolved model-predicted SMFs (see Figure 13), in the red- 
shift range 1.3 < z < 2.0, where the high-mass end is rea- 
sonably well reproduced by all but the model of Wang et al. 
(2008), all models significantly over-predict the number den- 
sity of galaxies below the characteristic stellar mass. This is 
true also in the redshift range 2.0 < z < 3.0, but now all mod- 
els also significantly under-predicting the number densities of 
the massive galaxies. The disagreement at the high-mass end 
is even more pronounced in the redshift range 3.0 < z < 4.0, 
where aU models show a significant deficiency of massive 
galaxies with respect to observations. If the comparison be- 
tween the observed and the model-predicted SMFs is instead 
performed using the convolved model-predicted SMFs, in- 
tended to represent measurement errors in logM^tar (see Fig- 
ure 14), the disagreements at the high-mass end are signifi- 
cantly reduced, although the model of Somerville et al. (2008) 
now tends to over-predict the number density of galaxies in 
the redshift range 1.3 < z < 2.0 at all stellar masses, while 
the model of Monaco et al. (2007) still significantly under- 
predicts the number density of massive galaxies at 3.0 < z < 
4.0. Convolving the model-predicted SMFs does not instead 
help at all to solve the discrepancies below the characteristic 
stellar mass. These results are in qualitatively agreement with 
the comparison of observed and predicted rest-frame optical 
luminosity functions performed at 2.0 < z < 3.3 by March- 
esini & van Dokkum (2007), who found that aU models signif- 
icantly over-predict the observed number density of galaxies 
at the faint-end. 

As already pointed out in Somerville et al. (2008), poten- 
tially serious discrepancies, common to all of the CDM-based 
semi-analytic models, are connected with low-mass galaxies. 
This is indeed what we find, although the presence of a sig- 
nificant population of very massive galaxies out to redshift 
2^3-4 and the little observed evolution in its number densi- 
ties from z = 4.0toz = 1-3 highUght other potential problems 
within the theoretical models. 

The different comparison between observed and predicted 
SMFs in the two stellar mass regimes (below and above 
Mstai- ^ 10" Mq) is also shown in Figures 15 and 16, where 
the evolution of Pstai ^^'^ pI^^'^^'^ function of redshift 
are compared to the predictions from the semi-analytic mod- 
els before and after convolution of the model-predicted SMFs 
with a normal distribution of standard deviation 0.25 dex, re- 
spectively. 

If no convolution with a normal distribution of standard 
deviation 0.25 dex is applied to the model-predicted SMFs 
(Figure 15), aU models severely under-predictthe stellar mass 
density of the massive galaxies at all redshifts. Only the 
model of Somerville el al. (2008) well match the stellar mass 
density of massive galaxies at 1.3 < z < 2.0, although it suf- 
fers from the same deficit of massive galaxies at z > 2.0. As 
in the comparison with the SMFs, if the convolved model- 
predicted SMFs are adopted (Figure 16), the disagreements 
for the stellar mass densities of massive galaxies is signifi- 
cantly reduced, especially for the model of Wang et al. (2008) 
which now reasonably well reproduces the evolution of the 
stellar mass density of massive galaxies over the entire red- 
shift range 0.0 < z < 4.0. However, the model of Monaco 
et al. (2007) still under-predicts significantly the stellar mass 
density of massive galaxies at 3.0 < z < 4.0, while the model 
of Somerville et al. (2008) now over-predicts the stellar mass 
density of massive galaxies at 1.3 < z < 2.0. The moderate 
success (after convolution) at the high-mass end is however 
counter-balanced by the failure in predicting the evolution of 
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11 10 11 10 11 10 11 10 

Log(M„„ [Ms„„]) Log(M„„ [M^„]) Log(M„„ [Ms„„]) Log(M„„ [Ms„„]) 

Fig . 13. — Comparison between the observed SMFs and the semi-analytic model predicted SMFs not convolved with a normal distribution of standard deviation 
0.25 dex. The predicted SMFs are represented by black curves: solid for the Somerville et al. (2008) model, dotted for the Monaco et al. (2007) model, and 
dashed for the Wang et al. (2008) model. Red, green, and blue filled circles represent the SMFs measured from the 1 /Vmax method; the thick eiTor bars include 
Poisson errors, cosmic variance, and the uncertainties from photometric redshift random errors; the thin eiTor bars include the systematic uncertainties, with the 
exclusion of the effect due to the use of the bottom-light IMFs for consistency with the theoretical models. Red, green, and blue solid curves represent the SMFs 
measured from the maximum-likelihood analysis. The shaded regions represent the total 1 cr en'ors, including the systematic uncertainties. The top panels show 
the comparisons between the observed and the model-predicted SMFs. The bottom panels show the differences between the model-predicted SMFs and those 
derived in this work, A$ = log "3? models - log 'I'ours, plotted as function of stellar mass. The SMFs predicted from the models are in general too steep, significantly 
over-predicting the number densities of galaxies below the characteristic stellar mass at all redshifts, and severely under-predicting the number density of the 
massive galaxies at z > 2.0. The good agreement at z = 0. 1 is the direct result of the optimization of the free parameters in the models to match the z = universe. 




11 10 11 10 11 10 11 10 

Log(M.,„ [Ms„„]) Log(M„„ [M^„]) Log(M„„ [Ms„„]) Log(M.,„ [M^„]) 

Fig. 14. — Same as in Fig. 13, except here the SMFs predicted by the semi-analytic models have been convolved with a normal distribution of standard 
deviation 0.25 dex; no convolution has been applied to the model-predict SMF at z = 0.1. While the discrepancies at the low-mass end are still present, the 
discrepancies at the high-mass end are significantly reduced. 



22 



Marchesini et al. 





8.5 






u 

Q. 


8.0 








7.5 




7.0 


a« 




o 

—I 


6.5 



1 1 1 1 


1 1 1 1 

10'°<M<10" 


1 1 1 

10"<M<10 


1 

12 . 
















V, X 

■ 

v., \* 


NOT CONVOLVED | 








V, 1 \ — f — ' 

\ T ' — 




SOB 








t \ 1 

x \ 




M07 








\ \ 




W&DL08 








.\ \. 




1 2 3 




1 2 3 


1 2 3 


z 




2 


z 



Fig. 15. — Comparison of the observed stellar mass densities with those 
predicted by the semi-analytic models as function of redshift, with no con- 
volution with a normal distribution of standard deviation 0.25 dex appUed to 
the model-predicted SMFs. Filled circles represent the observed stellar mass 
densities, while the curves are those predicted by the models. The left, mid- 
dle, and right panels show the comparison between observations and model 
predictions for pl;^^^\ p'S?'^'^". md plti?^^'^ respectively The thick 
error bars include Poisson errors, cosmic variance, and the uncertainties from 
photometric random errors; the thin error bars include the systematic un- 
certainties, with the exclusion of the effect due to the use of the bottom- 
Ught IMFs. The models fail in reproducing the observed evolution of the 
stellar mass density of low-mass galaxies (10'" < Msu^/Mq < 10"), and 
severely under-predict the stellar mass density of massive galaxies (10" < 
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Fig. 16. — Same as Fig. 15, except here the predictions of the semi- 
analytic models are convolved with a normal distribution of standard devi- 
ation 0.25 dex. While all models fail in reproducing the observed evolution 
of the stellar mass density of low-mass galaxies (10'° KM^i^/Mq < lO"), 
the disagreements between the observed and predicted stellar mass densities 
of massive galaxies are significantly reduced. 

the stellar mass density of galaxies below the characteristic 
stellar mass (10'° < Mstar/M© < lO"). Over the entire red- 
shift interval 1.3 < z < 4.0, all models over-predict the ob- 
served stellar mass density of low-mass galaxies by a factor of 
a few, with the predictions from the Somerville et al. (2008) 
model showing the largest disagreements with the observa- 
tions. 

We conclude that the models succeed for some redshift 
and stellar mass ranges and fail for others, although this may 
partly be due to systematic errors in the observed SMFs and 
stellar mass densities (see § 5 and § 6). In general, the 
model-predicted SMFs are too steep compared to the ob- 
served SMFs, resulting in a significant over-prediction of the 
number densities of low-mass galaxies. Models also fail in 
predicting the observed number density of massive galaxies 
observed at 3.0 < z < 4.0. 

8. SUMMARY AND CONCLUSIONS 

In this paper we have measured the stellar mass functions 
of galaxies at redshifts 1.3 < z < 4.0 from a composite K- 
selected sample constructed with the deep NIR MUSYC, the 
ultra-deep FIRES, and the GOODS-CDFS surveys, all having 
very high-quaUty optical-to-MIR data. This sample is unique 



in that it combines data from surveys with a large range of 
depths and areas in a self-consistent way. 

The total effective surveyed area of ^ 5 1 1 arcmin^ over sev- 
eral independent field of views allowed us to minimize the un- 
certainties due to cosmic variance, and empirically quantify 
its contribution to the total error budget. Moreover, we were 
able to probe the high-mass end of the SMFs with unprece- 
dented good statistics. The three adopted surveys allowed 
us to empirically derive the redshift-dependent completeness 
Umits of the ^-selected samples by exploiting their different 
depths. This is a significant improvement with respect to pre- 
vious studies, since it does not rely on stellar population syn- 
thesis models to assess the completeness in stellar mass of the 
sample. Finally, the ultra-deep FIRES survey allowed us to 
probe the low-mass end of the SMF down to stellar masses as 
small as ~ 0.05 times the characteristic stellar mass. 

We provide, for the first time, a comprehensive analysis of 
all random and systematic uncertainties affecting the derived 
SMFs. We have quantified the uncertainties on the SMF due 
to photometric redshift random errors by repeating the full 
analysis on a set of 100 mock A'-selected catalogs created by 
perturbing the observed photometry of each object according 
to its formal errors. Systematic uncertainties due to differ- 
ent sets of SED-modehng assumptions and different input pa- 
rameters in the estimate of the photometric redshifts have also 
been quantified by changing IMF, metalUcity, stellar popula- 
tion synthesis model, extinction curve in the SED-modeling, 
and template set and template error function in EAZY, the 
code used to estimate photometric redshifts. 

We found that the stellar mass density evolves by a fac- 
tor of ~ 17^jQ from z = 4.0, and a factor of ~ 4 ± 0.2 since 
z = 1.3. The observed evolution appears to be mostly driven 
by a change in the normalization $* of the SMF, rather than a 
change in the slope a or the characteristic stellar mass M*^^^, 
especially since redshift z = 3.0. This is by itself a very in- 
teresting result, since it implies that the physical processes re- 
sponsible in the building up of the stellar content of galaxies 
since z = 3.0 do not change significantly the shape of the SMF 
(defined by a and M*^^). The shape of the SMF is instead 
quite different at z ~ 3.5, with a much steeper low-mass end 
slope. We note however that the SMF at 3 .0 < z < 4.0 is quite 
uncertain. We also find evidence for mass-dependent evolu- 
tion of the SMF of galaxies with redshift. Specifically, galax- 
ies below the characteristic stellar mass show stronger evolu- 
tion with cosmic time with respect to the galaxies at the high- 
mass end, with the most massive galaxies (Mstar > lO''-^ Mq) 
showing a remarkable lack of evolution. 

If both random and systematic errors affecting the observed 
SMFs are taken into account, the previous results are no 
longer robust, and significant progress is required to better 
constrain the high-redshift SMFs. The observed absolute 
number of the most massive galaxies is very small, only a 
handful in each of the targeted redshift interval. Surveys with 
much larger area are needed in order to significantly improve 
the statistics on these rare objects. Cosmic variance, i.e., 
field-to-field variations, remains a significant source of un- 
certainty, especially at the high-mass end and at z < 2. There- 
fore, large surveys over multiple spatially disjoint fields are 
required to make significant progress in this respect. Last, but 
not least, most of the redshift information comes from pho- 
tometric redshift estimates. At the very high-mass end, pho- 
tometric redshift errors are a non-neghgible source of errors, 
since the uncertainties wiU tend to scatter galaxies from the 
well-populated region around M*^^ toward the sparsely popu- 
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lated high-mass end. Redshift errors are a significant source 
of uncertainties also at the highest considered redshift inter- 
vals. 

Obtaining large numbers of spectroscopic redshifts for K- 
selected high-z sources has proven difficult and extremely 
time consuming, even for very bright galaxies (Kriek et al. 
2008). Note also that the galaxies belonging to the high-mass 
end (Mstar > 10^' M©) are generally very faint in the observed 
optical {R ^ 26; van Dokkum et al. 2006), hence making op- 
tical spectroscopy completely unfeasible. Although the situa- 
tion will be improved in the near future by the advent of multi- 
object NIR spectrographs, which will allow to construct large 
sample of high-z /T-selected galaxies with spectroscopic red- 
shift measurements, these projects will only probe the bright- 
est population, leaving the high-mass end of the SMF stiU 
prone to photometric redshift uncertainties affecting the rest 
of the galaxy population. 

The on-going NEWFIRM Medium-Band Survey (van 
Dokkum et al. 2009) will be able, upon completion, to im- 
prove on all of the above aspects. This deep and wide- 
field NOAO/Yale survey uses the newly commissioned NEW- 
FIRM instrument (mounted on the 4 m Mayall telescope) 
with custom-made medium band-width filters over the wave- 
length range 1-1.8 //m to obtain well-sampled SEDs and 
high quahty photometric redshifts for K^^^ < 23.4 galaxies. 
The full survey will provide accurate redshift measurements 
(Az/(1 +z) ~ 0.02) for - 8000 galaxies at 1 .5 < z < 3 .5 over a 
total area of --0.5 square degree in the COSMOS and AEGIS 
fields, allowing to drastically reduce the impact of random 
uncertainties and cosmic variance in the measurements of the 
high-redshiftSMFs. 

We also note that we cannot constrain the low-mass end of 
the SMFs very well at z > 2. If the SMF is much steeper at 
the low-mass end than indicated by our maximum-likelihood 
analysis, we could be missing a significant fraction of the in- 
tegrated stellar mass (see Reddy & Steidel 2009). Reddy & 
Steidel (2009) suggest that up to --50% of the total stellar 
mass at 1.9 < z < 3.4 is in faint-UV galaxies with masses 
smaller than ^10"^ Mq (compared to ^10%-20% from an 
extrapolation of our Schechter fits). However, Reddy & Stei- 
del (2009) do not measure stellar masses directly, but convert 
UV luminosity to stellar mass. Ultra-deep near-IR imaging 
can directly determine whether there is a large population of 
very low-mass galaxies at high redshift which contribute sig- 
nificantly to the total stellar mass density. Deep surveys with 
WFC3 on HST, as well as deep ground-based surveys such as 
the UKTOSS Ultra-Deep Survey (Lawrence et al. 2007) and 
the Ultra- VISTA survey'^, will allow us to much better con- 
strain the low-mass end of the SMFs. 

Progress in observations have to be supported by signifi- 
cant progress in the theoretical arena. Systematic uncertain- 
ties due to different SED-modeUng assumptions are a signif- 
icant, if not dominant, contribution to the total error budget. 
Convergence on the stellar population synthesis models is of 
paramount importance not only to significantly decrease the 
systematic uncertainties on the derived SMF, but also to have 
a correct understanding of the properties of galaxies, such as 
age, star formation rate, etc. Understanding the IMF and its 
evolution as function of cosmic time is also extremely impor- 
tant, especially for studies at high-redshift. Significant work 
is however required theoretically to better understand what 



physical parameter is most responsible in shaping the IMF 
and its evolution with cosmic time. 

The observed SMFs have been compared with the SMFs 
predicted from the semi-analytic models of Monaco et al. 
(2007), Wang et al. (2008), and Somerville et al. (2008). The 
model-predicted SMFs are generally too steep with respect to 
the observed SMFs, resulting in significant over-prediction of 
the number densities of galaxies at the low mass end. While 
relatively good agreement is observed at the high-mass end 
at z < 2.5, some models tend to under-predict the observed 
number densities of massive galaxies at 2.5 < z < 4.0. The 
discrepancy at the high-mass end is susceptible to uncertain- 
ties in the models and the data, but the discrepancy at the low- 
mass end may be more difficult to explain. These results are 
robust, even when all random and systematic uncertainties are 
included, suggesting that current models do not yet provide a 
complete description of galaxy formation and evolution since 
z = 4.0. 
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TABLE A 1 

Characteristics of the IRAC observations 



rilter 


h-xposure Time 


r WriJVl 


Total Limiting Magnituue 


Positional Accuracy 


Oalactic extinction 


(jim) 








V_cUCaCC J 


(mag) 


SDSS-1030: 












3.6 


1.0 


1.4 


24.53 


0.07 


0.005 


4.5 


1.0 


1.4 


24.14 


0.08 


0.004 


5.8 


1.0 


1.7 


22.35 


0.07 


0.004 


8.0 


1.0 


2.0 


22.16 


0.08 


0.003 


CW-1255: 












3.6 


1.0 


1.5 


24.74 


0.05 


0.003 


4.5 


1.0 


1.4 


24.37 


0.07 


0.003 


5.8 


1.0 


1.8 


22.52 


0.05 


0.002 


8.0 


1.0 


2.1 


22.45 


0.07 


0.002 


HDFS-2: 












3.6 


0.5 


1.4 


24.49 


0.08 


0.005 


4.5 


0.5 


1.4 


24.00 


0.09 


0.005 


5.8 


0.5 


1.8 


22.21 


0.08 


0.004 


8.0 


0.5 


2.2 


22.08 


0.09 


0.004 


HDFS-1: 












3.6 


0.3-4.2 


1.7 


24.41-25.06<^ 


0.06 


0.005 


4.5 


0.3-4.2 


1.7 


24.12-24.83<^ 


0.07 


0.005 


5.8 


0.3-4.2 


2.2 


22.39-23.36"^ 


0.06 


0.004 


8.0 


0.3-4.2 


2.1 


22.35-23.24<^ 


0.07 


0.004 



"The total limiting magnitudes were estimated using the empty aperture method, corrected for the flux 
missed outside the 3" aperture used for photometry. This correction is ~ 0. 18 mag. 

''The rms difference between bright star positions in IRAC and A'-band image, after pointing refinement. 

"^The first number corresponds to the shallower area, while the second depth corresponds to the deeper area 
overlapping with the HDF-S Proper field. 

APPENDIX 

DATA REDUCTION AND PHOTOMETRY OF THE SPITZER-IRAC DATA IN MUSYC: 

Here we describe the data taken with the Infrared Array Camera (IRAC; Fazio et al. 2004) on board Spitzer over the deep NIR 
MUSYC fields, their reduction, and the creation of /T-selected catalogs with IRAC photometry. The /T-selected catalogs with 
IRAC photometry included is publicly available at http://www.astro.yale.edu/musyc. 

Spitzer-IRAC data 

The IRAC data over the deep NIR MUSYC fields come from two different sources. Specifically, the IRAC data over the 
HDFS-2, the SDSS-1030, and the CW-1255 fields come from the Spitzer Space Telescope Cycle-3 program GO-30873 (RI.: 
Labbe), while the IRAC data over the HDFS-1 field are part of the GTO-214 program (RL: Fazio). Table Al summarizes the 
characteristics of the IRAC data over the deep MUSYC fields, such as the total exposure time, full-width half maximum (FWHM), 
limiting depth, positional accuracy, and Galactic extinction in each of the four IRAC bands. The total exposure times vary from 
30 min for the HDFS-2 field, to 1 hr in the SDSS-1030 and CW-1255 fields across the entire 10' x 10' field. The exposure time 
across the HDFS-I is instead not homogeneous, being ~ 4.2 hrs in the small area overlapping with the HDF-S Proper field, and 
20 min everywhere else. All fields were covered in all four IRAC channel, namely the 3.6, 4.5, 5.8, and 8.0 fim bands. 

Data reduction 

The reduction started with the basic calibrated data (BCD) as provided by the Spitzer Science Center pipeline. We applied a 
series of procedures to reject cosmic rays and remove artifacts such as column pulldown, muxbleed, maxstripe, and the "first- 
frame effect" (Hora et al. 2004). Then, the single background-subtracted frames were combined into a mosaic image large enough 
to hold all input frames (using WCS astrometry to align the images). In this step, bad pixels were also masked, and distortion 
corrections applied in WCS. This image, in combination with the /T-band reference image, was then used to refine the pointing of 
the individual mosaics (6 mosaics, in a 2 x 3 grid). After pointing refinement, positional accuracy is ~ 0.05-0.09 arcsec. We note 
that the source-fitting algorithm developed by I. Labbe et al. (in preparation) that we have used to derive the IRAC photometry 
takes care of residual shifts. Finally, the individual pointing-refined frames were registered to and projected on the public /T-band 
images (0.27" pixel scale) of the four MUSYC fields (Quadri et al. 2007),'' and average-combined. Flux conservation has been 
forced throughout the reduction. 

Seeing, zero-points, and limiting depths 

FWHM of the IRAC images were derived by taking the median of the FWHM for a set of ~ 15-30 bright, isolated stars in the 
fields. The FWHM of each field and each IRAC band are listed in Table Al, and amount to 1.4-1.7", 1.4-1.7", 1.7-2.2", and 
2.0-2.2" for the 3.6, 4.5, 5.8, and 8.0 /xm bands, respectively. 

AH optical and NIR data and the *r-selected catalogs from the deep NIR MUSYC survey are publicly available from http://www.astro.yale.edu/musyc. 
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Following Quadri et al. (2007), the IRAC photometry in the J?'-selected catalogs are presented in units of flux, normalized 
so that the zero-points is 25 on the AB system. The use of flux, rather than magnitudes, avoids the problem of converting the 
measured flux uncertainties into magnitude uncertainties, the problem of asymmetric magnitude uncertainties for low S/N objects, 
and the loss of information for objects that have negative measured fluxes. Vega to AB magnitude conversions are 2.78, 3.26, 
3.75, and 4.38 for the 3.6, 4.5, 5.8, and 8.0 /xm bands, respectively. We adopted the following transformation from Jy to Vega 
magnitudes: 

wvega = -2.5 log iflux[Jy]) +B, (Al) 

with B = 6.12, 5.64, 5.15, and 4.52 for the 3.6, 4.5, 5.8, and 8.0 fim bands, respectively. The IRAC zero-point uncertainties are 
of the order of 2% (Reach et al. 2005), which was included in the flux error budget. 

The IRAC photometry has been corrected for Galactic extinction, using the extinctions listed in Table Al. These values, 
taken from the Galactic Dust Extinction Service (http://irsa.ipac.caltech.edu/applications/DUST), were derived using the data 
and technique of Schlegel et al. (1998). 

The properties of the IRAC images were analyzed following the same approach as for the deep NIR MUSYC data (Quadri et 
al. 2007). Briefly, the technique uses aperture photometry distributed randomly over empty regions of the image to quantify the 
rms of background pixels within the considered aperture size. For a given aperture size, the distribution of empty aperture fluxes 
is well fitted by a Gaussian. The total limiting AB magnitudes (3 a for point sources, i.e., corrected for the flux missed outside 
the aperture used for photometry) are Usted in Table Al . 

Photometry 

A source-fitting algorithm developed by I. Labbe et al. (2008, in preparation), especially suited for heavily confused images 
for which a higher resolution prior (in this case the /T-band image) is available, was used to extract the photometry from the 
IRAC images. A short description with illustration was also presented by Wuyts et al. (2007). Since this program does not take 
into account large-scale background variations, these were removed a priori from the IRAC images. Briefly, the information on 
position and extent of sources bases on the higher resolution /T-band segmentation map was used to model the lower resolution 
IRAC 3.6-8.0 /im images. Each source was extracted separately from the ^T-band image and, under the assumption of negligible 
morphological A:-corrections, convolved to the IRAC resolution using the local kernel coefficients. Convolution kernels were 
constructed using bright, isolated, unsaturated sources in the K and the IRAC bands (derived by fitting a series of Gaussian- 
weighted Hermite functions to the Fourier transform of the sources, rejecting outlying or poorly-fitting kernels), and a smoothed 
2D map of the kernel coefficients was stored. A fit to the IRAC image was then made for all sources simultaneously, where the 
fluxes of the objects were left as free parameters. Next, we subtracted the modeled light of neighboring objects and measured 
the flux on the cleaned IRAC map within a fixed 3" diameter aperture. Through a visual inspection of the IRAC residual 
image with all sources subtracted, we conclude that this method effectively removes contaminating sources (for an illustration 
of this technique, see also Figure 1 in Wuyts et al. 2007). In order to compute a consistent K-IRAC color, we measured the 
source's flux /com, k on a cleaned /T-band image convolved to the IRAC resolution within the same aperture. We then scaled the 
photometry to the same color apertures that were used for the NIR photometry, allowing a straightforward computation of colors 
over U-Xo-% jim wavelength baseline. For the IRAC photometry, this means that the catalog flux was computed as follows: 

/iRAC.col = /lRAC,3" J^^''^°^ . (A2) 
/conv,K,3" 

Note that the used source-fitting algorithm developed by 1. Labbe et al. (in prep.) takes into account the spatial extent of the 
sources on the reference K-hand image, and it does not adopt the best-fit flux from the source fitting as final photometry, but 
rather measures the flux within an aperture on the cleaned image (followed by an aperture correction). This allows more robust 
photometry in cases where the object profile varies from the reference to the low-resolution image. Uncertainties in the measured 
fluxes in the IRAC bands were derived as described in Wuyts et al. (2008), accounting for the background rms and residual 
contamination of the subtracted neighbors (I. Labbe et al. 2008, in prep.). 

ILLUSTRATION OF COMPLETENESS ESTIMATION: THE SDSS-1030 SAMPLE 

We have used a different approach (described in § 4. 1) to estimate the redshift-dependent completeness limit in stellar mass of 
our /^-selected sample to be used to derive the SMFs of galaxies. In the following, we illustrate the estimation of the redshift- 
dependent completeness limit in stellar mass for the SDSS-1030 sample. 

First, we selected galaxies belonging to the available deeper samples, namely the HDFS, the MS-1054, and the CDFS samples. 
Second, we scale their fluxes and stellar masses to match the SDSS-1030 /T-band 90% completeness limit. This is illustrated 
in Figure Bl. In this figure, the colored filled symbols represent objects from the deeper samples (HDFS in blue, MS-1054 in 
red, and CDFS in orange) scaled up in flux to the SDSS-1030 /T-band 90% completeness hmit. These objects represent objects 
immediately at our detection Umit. The upper envelope of these points in the (Ms,aj scaled -z) space represents the most massive 
galaxies that might escape detection/selection in our analysis. Therefore, this upper envelope, encompassing 95% of the points, 
provides a redshift-dependent stellar mass completeness limit for the SDSS-1030 sample. 

This suggests that at z ^ 2.3, the SDSS-1030 sample is approximately complete for stellar masses Mstar> 1O"M0. In the same 
panel, we also compare the empirically-derived completeness limit (solid curve) with the completeness derived assuming an SSP 
with no dust formed at z = 10 and scaled to match the SDSS-1030 ^-band 90% completeness limit (dashed curve). It is obvious 
from this comparison that, for SDSS-1030, the SSP-derived completeness is similar to the empirically-derived completeness only 
for z < 2, while at z > 2 the SSP-derived completeness implies a higher completeness in stellar mass than empirically derived. The 
difference between the empiricaUy- and the SSP-derived completeness is a function of redshift, with the differences increasing 



26 



Marchesini et al. 




1.5 2.0 2.5 3.0 3.5 4.0 
z 

Fig. B 1 . — Empirically-determined completeness in stellar mass as a function of redshift; black filled circles represent the stellar masses for the SDSS-1030 
galaxies scaled down in flux to match the SDSS-1030 A'-band 90% completeness limit, and plotted as a function of redshift. The other circles show the stellar 
masses for the CDFS (orange), the MS-1054 (red), and the HDFS (blue) galaxies scaled in flux to match the SDSS-1030 S'-band 90% completeness limit. 
The solid curve represents the upper envelope of these points (encompassing 95% of the points), and effectively defines, as a function of redshift, the limiting 
stellar mass corresponding to the observed flux limit of the SDSS-1030 sample. At z ^ 2.3, the SDSS-1030 sample is approximately complete for stellar masses 
Mstai- > lO" . We emphasize that the stellar masses plotted here are not the actual stellar masses, but the stellar masses scaled to the A'-band limit of the SDSS-1030 
field. 

going to higher redshift. This could be due, i.e, to the fact that dust extinction becomes progressively more important for galaxies 
in the high-mass end with increasing redshift. 

COMPARISON WITH PREVIOUSLY PUBLISHED GALAXY STELLAR MASS FUNCTIONS: 

Here we compare our results to previous studies on the SMFs of galaxies at z > L These works include the work of Drory 
et al. (2005), Fontana et al. (2006), Pozzetti et al. (2007), Eisner et al. (2008), and Perez-Gonzalez et al. (2008). None of these 
works have included a comprehensive analysis of the uncertainties (random and systematic) on the derived SMFs. As already 
pointed out, our work represents the first analysis of the SMFs at high-redshift with a comprehensive analysis of the errors. The 
statistically significant disagreements among the different works mostly stem from the lack of a complete analysis of the errors 
of the SMFs from the literature . Once a complete analysis of the errors is performed, as done for the first time in our work, the 
disagreements between the different measurements of the SMFs are no longer statistically significant. 

Our /T-selected composite sample is unique in that it combines very high-quality multi-waveband data from surveys with a large 
range of depths and areas in a self-consistent way. The available large number of single large area fields allows for the empirically 
estimate of the contribution of cosmic variance to the error budget, which is the dominant source of random errors in the lowest 
targeted redshift range and at the high-mass end. The combination of deep and ultra-deep samples allows us to empirically derive 
the completeness in stellar mass of the ^T-selected sample, without relying on stellar models, such as a passively evolving single 
stellar population formed at very high redshift (z ~ 10) without any dust extinction (SSP-derived completeness), as done in most 
of the studies in the literature. 

We therefore conclude that, with respect to the works in the literature, our SMFs are a significant improvement due to 1) the 
comprehensive analysis of the errors of the SMFs, both random and systematic (completely missing in the literature); 2) the large 
range in stellar masses probed by the composite /T-selected sample, which allowed for better sampling of both the high- and the 
low-mass end; 3) the large surveyed area and the large number of independent fields, allowing for a decrease of the uncertainty 
due to cosmic variance (with respect to previous works) and to empirically quantify its contribution; 4) the availability of samples 
with different depths, allowing us to empirically derive, for the first time, the redshift-dependent completeness limits in stellar 
mass. 

Drory etal (2005): 



Drory et al. (2005) derived the SMFs of galaxies from z = 0.25 to z = 5 from two different samples: an /-selected sample over 
40 arcmin^ of the FORS Deep Field (FDF; Heidt et al. 2003) with UBgRIzJK coverage, consisting of 5557 galaxies down to 
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Fig. CI.— Comparison between the SMFs from Drory et al. (2005) and from this work at 3.0 < z < 4.0 (left), 2.0 < z < 3.0 (middle), and 1.3 < z < 2.0 
(right). Symbols as in Fig. 11, except black filled symbols and black curves representing the SMFs derived in Drory et al. (2005). The SMFs of Drory et al. 
(2005) have been scaled by -0.2 dex along the x-axis to take into account difference in the adopted IMF. Note that the SMFs from Drory et al. (2005) are actually 
derived in the following redshift intervals: 3.00 < z < 4.00 (left panels), 2.25 < z < 3.00 (middle panels), 1.25 < z < 1.75 and 1.75 < z < 2.25 (right panels; 
filled black circles and squares, respectively). 

/ 26.8 (50% completeness); and a /T-selected sample over 50 arcmin^ in the GOODS-CDFS field with UBVRIJHK coverage. 
No mid-IR IRAC photometry was used. The quoted photometric redshift accuracy in Az/{1 + Zspec) is ^ 0.03. The stellar 
masses were derived by fitting the observed SEDs with a model grid of BC03 models. The star-formation history (SFH) was 
parameterized by a two-component model, with a main component with a smooth SFH modulated by a burst of star formation. 
The main component is parameterized by an exponentially declining SFH with an e-folding timescale r G [0.1, oo] Gyr, and a 
metallicity of -0.6 < [Fe/H] < 0.3. The age was allowed to vary between 0.5 Gyr and the age of the universe at the object's 
redshift. This component was linearly combined with a burst modeled as a 100 Myr old CSF rate episode of solar metallicity, 
restricting the burst fraction Q < (3 < 0.15 in mass. A Salpeter (1955) IMF truncated at 0.1 and 100 M0 was adopted for both 
components. Finally, both components are allowed to exhibit a different and variable amount of extinction by dust. SMFs were 
derived with the 1 /Vmax method, although it is unclear how the completeness is stellar mass as a function of redshift was derived. 
The Schechter function parameters were derived by fitting the 1 /Vmax points with a Schechter function. 

Their surveyed area is sim6.5 times smaller than the surveyed area of our work, making their derived SMFs very much affected 
by cosmic variance. A direct comparison between the SMFs of Drory et al. (2005) and the SMFs measured in our work is shown 
in Figure CI. Note that the redshifts bins used in Drory et al. (2005) are not exactly the same as used in the present work. 

The top panels of Figure CI show the SMFs derived in our work without the inclusion of cosmic variance and systematic 
uncertainties in the plotted error bars; error bars include Poisson errors and photometric redshift random errors. The errors of 
the Drory et al. (2005) SMFs only include random uncertainties. From the top panels of Figure CI, the SMFs from Drory et al. 
(2005) are consistent with those derived in our work only at the high-mass end, while their number densities are significantly 
higher at the low-mass end. The slopes of the SMFs at the low-mass end are steeper than those derived in our work, especially 
at z '--^ 1 .6. In the highest redshift range, where Poisson statistics dominates the error budget, the two SMFs are consistent within 
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2 (7 in the overlapping stellar mass regime. 

The SMFs of Drory et al. (2005) become fully consistent with ours at the high-mass end once all sources of errors are included 
in the error budget (bottom panels of Figure CI). However, the discrepancy at the low-mass end is still present and significant. 
The much higher densities derived by Drory et al. (2005) at the low-mass end (hence, the steeper slope a) are very likely caused 
by the use of the 1 /Vmax method in combination with an inappropriate redshift-dependent completeness in stellar mass. While 
it is unclear how the completeness in stellar mass has been derived in Drory et al. (2005), we notice that, if an SSP-derived 
completeness were to be adopted in place of the empirically-derived completeness, steeper densities at the low-mass end would 
be derived, due to the more conservative completeness of the former with respect to the latter. We therefore believe that the SMFs 
of Drory et al. (2005) at the low-mass end are systematically too large, perhaps due to the use of incorrect redshift-dependent 
completeness limits in stellar mass. 

Fontanaetal. (2006): 

Fontana et al. (2006) derived the SMFs of galaxies at 0.4 < z < 4.0 from the GOODS-MUSIC sample (Grazian et al. 2006). 
Their final /iT-selected sample consists of ^2931 galaxies (1762 with //-band coverage) complete down to Kg « 23.5, over an 
area of ^143.2 arcmin^ with UB4j,sV(,Q(,hi5Zii5{)JHK and IRAC bands coverage. Note that this dataset make use of the same public 
data of the FIREWORKS-CDFS dataset (Wuyts et al. 2008) that we have used in our work, but without the inclusion of the WFI 
BVRI band data. Consequently, the SEDs of sources in the FIREWORKS-CDFS catalog are better sampled. We also note that 
(Wuyts et al. 2008) pointed out a systematic offset in the IRAC fluxes of the GOODS-MUSIC catalog with respect to the IRAC 
fluxes of the FIREWORKS-CDFS catalog, with the former averagely fainter than the latter by ^ 30%. The observed offset in the 
IRAC photometry was largely attributed to the use of an early version of the IRAC PSF by Grazian et al. (2006) and a bug in the 
normalization of the smoothing kernel for the IRAC data by Grazian et al. (2006) (see Wuyts et al. 2008 for details). 

Stellar masses were derived by fitting the observed SEDs with a set of templates computed with the BC03 spectral synthesis 
models. A Salpeter (1955) IMF was adopted with various metalUcities (from Z = 0.02 Zq to Z = 2.5 Zq), dust extinction 
(0 < E(B-V) < 1.1) with a Calzetti et al. (2000) extinction curve, and e-folding timescales (r G [0.1, 15] Gyr) in an exponentially 
declining SFH. 

SMFs were derived with the standard 1 /Vmax formalism and the maximum-likelihood analysis assuming a Schechter function 
with Schechter parameters a, M*^^, and $* evolving with redshift (for a total of seven free parameters, three of which constrained 
by the local SMFs derived in Cole et al. 2001). A treatment has been included to correct for the incompleteness in mass at the 
faintest levels (see Fontana et al. 2004 for details). Briefly, they start from the threshold computed from a passively evolving 
system (SSP-derived completeness limit), below which only a fraction of objects of given mass will be observed. Then, at 
any redshift, the observed distribution of M/L was obtained for objects close to the magnitude limit of the sample. Using this 
distribution, the fraction of lost galaxies as function of redshift and mass was computed. The correction is finally applied to 
the volume element Vmax of any galaxy in the 1 /Vmax biimed SMFs as well as in the number of detected galaxies entering the 
maximum-likelihood analysis. We note that this method assumes that the sample is complete for objects with stellar masses 
larger than the SSP-derived limit. This assumption is not necessarily valid at all redshifts, and depends on the specific depth of 
the sample (see the right panel of Figure 3). Moreover, their correction for incompleteness also assumes that the distribution of 
M/L ratio is independent of L. As shown in the left panel of Figure 3, this assumption might not be valid at all redshift, depending 
on the specific depth of the sample. In fact, at faint luminosities, the distribution of M/L is visibly different from that at the bright 
end. While it is hard to predict the differences on the low-mass end of the SMFs derived with their method with respect to 
the method used in our work (empirically-derived completeness limit), we notice that their correction has been applied only for 
sources for which the correction factor is smaller that 0.5 (usually affecting only the last 1 /Vmax point). Moreover, as shown in 
the right panel of Figure 3, the SSP-derived completeness Umit of the CDFS sample is quite similar to the empirically-derived 
completeness Umit. Therefore, we do not expect large differences in the derived SMFs at the low-mass end, except for the lowest 
stellar mass bins. 

A direct comparison between the SMFs of Fontana et al. (2006) and the SMFs measured in this work is shown in Figure C2. 
Again, top panels do not show the contribution of cosmic variance and systematic uncertainties in our SMFs, while they are 
included in the bottom panels of Figure C2. Error bars in the SMFs of Fontana et al. (2006) include Poisson errors and photometric 
redshift uncertainties, but no cosmic variance nor systematic uncertainties due to different SED-modeling assumptions. Note that 
the area surveyed in the GOODS-MUSIC sample is a factor of ~ 4 smaller than the surveyed area in our work, and it consists of 
a single pointing, making their SMFs significantly affected by cosmic variance. 

First of, as shown in the top panels of Fig. C2, the SMFs estimated from the 1 /Vmax method are broadly consistent within the 
errors. There appears to be a systematic offset, such that the SMFs of Fontana et al. (2006) are systematically shifted to lower 
masses. Although it is not simple to fully understand the origin of this difference, it could partly be due to the found systematic 
offset in the IRAC photometry. The SMFs derived using the maximum-likelihood analysis show a larger degree of disagreement, 
especially in the high-mass end at z 3.5. We stress once again that cosmic variance is a dominant source of uncertainty in the 
SMFs of Fontana et al. (2006) since they are derived from a single, relatively small field. The found differences at the bright end 
can well be accounted for by field-to-field variations. 

From the bottom panels of Figure C2, it is obvious that the SMFs of Fontana et al. (2006) are fully consistent with ours once all 
the source of uncertainties are taken into account, including the systematic effects due to different SED-modeling assumptions. 

Pozzetti et al. (2007): 

Pozzetti et al. (2007) derived the SMFs of galaxies from z = 0.05 to z = 2.5 from the ^T-selected sample of the VIMOS-VLT 
Deep Survey (VVDS; Le Fevre et al. 2005) 02h field. Their X'-selected sample consists of a shallow and a deeper component. 
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Fig. C2.— Comparison between the SMFs from Fontana et al. (2006) and from this worlc at 3.0 <z< 4.0 (left), 2.0 < z < 3.0 (middle), and 1.3 < z < 2.0 
(right). Symbols as in Fig. CI. The SMFs of Fontana et al. (2006) have been scaled by -0.2 dex along the jc-axis to take into account difference in the adopted 
IMF. Note that the SMFs from Fontana et al. (2006) are actually derived in the following redshift intervals: 3.0 < z < 4.0 (left panels), 2.0 < z < 3.0 (middle 
panels), 1.3 < z < 1.6 and 1.6 < z < 2.0 (right panels; filled black circles and squares, respectively). 

The shallow component consists of 6720 galaxies at < z < 2.5 down to K = 22.34 (90% complete) over 442 arcmin^, while the 
deeper component is made of 3440 galaxies down to K = 22.84 (90% complete) over 172 arcmin^. About 15% of the galaxies 
have secure spectroscopic identification. The waveband coverage consists of UBVRIugrizJK, but no IRAC coverage (note that 
only 170 arcmin^ have deep J and K coverage; lovino et al. 2005). Therefore, their surveyed area is ^^3 times smaller than the 
area surveyed in our work down to the same /T-band limit of K = 22.8. 

The stellar masses were derived by fitting the observed SEDs with a model grid of BC03 models. The star-formation history 
(SFH) was parameterized with an exponentially declining SFH with an e-folding timescale t e [0.1, oo] Gyr, and solar metallicity. 
The age was allowed to vary between 0.1 Gyr and the age of the universe at the object's redshift. A Chabrier (2003) IMF was 
adopted, and extinction modeled using the Calzetti et al. (2000) curve with Ay G [0,2.4]. 

SMFs were derived with the standard l/Vmax formalism and the maximum-likelihood analysis assuming a Schechter (1976) 
function. While Pozzetti et al. (2007) stress the important of using only galaxies with stellar masses above the stellar mass limit 
where all the SEDs are potentially observable (very restrictive limit), they ended up using as a lower limit of the mass range 
the minimum mass above which late-type SEDs are potentially observable. Therefore, we caution about potential biases in the 
estimate of the low-mass end of their SMFs. 

A direct comparison between the SMFs of Pozzetti et al. (2007) and the SMFs measured in our work is shown in Figure C3. 
Note that the redshift bins used in Pozzetti et al. (2007) are not exactly the same as used in the present work. The SMFs derived 
in our work are shown in the left panel without the contribution of cosmic variance and the systematic uncertainties, and in the 
right panel with all errors included. 

As shown in the left panel of Figure C3, the SMF from Pozzetti et al. (2007) at 1 .6 < z < 2.5 is perfectly consistent with ours. 
On the contrary, their SMF at 1.2 < z < 1.6 is significantly different at the high-mass end. If the two SMFs at 1.2 < z < 1.6 and 
1.6 < z < 2.5 are averaged, the obtained SMF is in very good agreement with the SMF derived in our work, except at the very 
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Fig. C3. — Comparison between the SMFs from Pozzetti et al. (2007) and from this work at 1.3 < z < 2.0. Symbols as in Fig. CI. Note that the SMFs from 
Pozzetti et al. (2007) are actually derived at 1.2 < z < 1.6 (tilled black squares) and 1.6 < z < 2.5 (filled black circles). 

low-mass end. We note that the low-mass end becomes consistent within the errors when the systematic uncertainties are taken 
into account in the total error budget, highlighting once again the importance of a comprehensive analysis of the errors. Overall, 
the agreement between the SMF of Pozzetti et al. (2007) and ours is very good. 

Eisner etal. (2008): 

Eisner et al. (2008) have derived SMFs of galaxies from z = 5.0 to z = 0.25 from the GOODS-MUSIC z + ZT-selected catalog 
of Grazian et al. (2006), the same used in Fontana et al. (2006). The catalog used in Eisner et al. (2008) is to first order zgso- 
band selected sample, with the addition of the remaining TiT-band sources that are undetected in the zgso-band, making it less 
trivial to understand the completeness of the sample, especially in stellar mass. This sample comprises 14800 galaxies down to 
Ziim ~ 26.0 (90% completeness level) over a total area of 143.2 arcmin^. 

The stellar masses were derived by fitting the observed SEDs with a model grid of BC03 models. The star-formation history 
(SFH) was parameterized with an exponentially declining SFH with an e-folding timescale t £ [0.5, 20] Gyr, and solar metallicity. 
The age was allowed to vary between 0.2 Gyr and the age of the universe at the object's redshift. A Salpeter (1955) IMF truncated 
at 0.1 and 100 Mq was adopted. Dust extinction was modeled using the Calzetti et al. (2000) curve with Ay €E [0, 1.5]. In addition 
to this main component, a starburst was superimposed which was allowed to contributed at most 20% to the z-band luminosity in 
the rest-frame. This component was modeled as a 50 Myr old episode of constant star formation with an independent extinction 
up to Ay = 2.0 mag. 

SMFs were derived with the 1 /Vmax method after correcting the data points for incompleteness due to the flux-limited sample. 
The completeness limit in stellar mass was estimated by scaling the z-band completeness limit by the calculated 95% quantile 
in M/Lz as a function of redshift, e.g. the limit below which 95% of the M/L^ ratios of the sample are located. The Schechter 
function parameters were then derived by fitting the 1 /Vmax points with a Schechter function after fixing the low-mass end slope 
at its error-weighted mean value. 

A direct comparison between the SMFs of Eisner et al. (2008) and the SMFs measured in our work is shown in Figure C4. 
Note that the redshifts bins used in Eisner et al. (2008) are not exactly the same as used in the present work. 

The SMFs estimated from Eisner et al. (2008) are in general good agreement with the SMFs measured in our work, especially 
at z 3.5 and z ~ 2.5. At the lower targeted redshift interval (1.3 < z < 2.0), the SMF of Eisner et al. (2008) shows a higher 
number density for the most massive galaxies with respect to our measurements. We note however that the sample of Eisner et 
al. (2008) is constructed from the single, relatively small GOODS field (^4 times smaller than the area surveyed in our work). 
Therefore, the SMFs derived by Eisner et al. (2008) are significantly affected by field-to-field variations, especially at low redshift 
and at the high-mass end. Note that the error bars in the SMFs of Eisner et al. (2008) do not include the error due to cosmic 
variance, which we have shown being the dominant contribution to the total random error budget at z < 2. We therefore conclude 
that the SMFs of Eisner et al. (2008) are fully consistent with our measurements, and that the disagreements at the high-mass end 
can be fully accounted by field-to-field variations. 

We finally note the large significant discrepancies between the SMFs derived from Eisner et al. (2008) and from Fontana et 
al. (2006). The used catalog is exactly the same, i.e., the GOODS-MUSIC catalog. The differences are in the way the stellar 
mass completeness limits are estimated (affecting the low-mass end of the SMF) and in the assumptions of the SED-modeling. 
Specifically, the metallicity in Eisner et al. (2008) is fixed to solar, while it is left as a free parameter in Fontana et al. (2006). 
Also, larger Ay are allowed in Fontana et al. (2006), as well as no secondary starburst component added, contrary to what done in 
Eisner et al. (2008). The SMFs from Fontana et al. (2006) are systematically lower at similar stellar mass, particularly at z ^ 3.5 
and at the high-mass end. While the differences are statistically significant, as already pointed out by Eisner et al. (2008), they 
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Fig. C4.— Comparison between the SMFs from Eisner et al. (2008) and from this work at 3.0 < z < 4.0 (left), 2.0 <z< 3.0 (middle), and 1.3 < z < 2.0 
(right). Symbols as in Fig. CI. The SMFs of Eisner et al. (2008) have been scaled by -0.2 dex along the x-axis to take into account difference in the adopted 
IMF. Note that the SMFs from Eisner et al. (2008) are actually derived at 3.01 < z < 4.01 (left panels), 2.25 < z < 3.01 (middle panels), 1.25 < z < 1.75 and 
1.75 < z < 2.25 (right panels; filled black circles and squares, respectively). 

can be due to the different choices of SED-modeling assumptions, stressing the importance of a comprehensive analysis of the 
errors, both random and systematics. 

Perez-Gonzalez et al. (2008): 

Perez-Gonzalez et al. (2008) derived the SMFs of galaxies from z = to z = 4 from a 3.6 /im and 4.5 /im Spitzer-IRAC selected 
sample. Their sample consists of ^19400 sources down to the 75% completeness limit (^--^23. 3 mag at 3.6 fj,m) over three fields, 
namely the HDF-N, the CDF-S, and the Lockman Hole fields, for a total surveyed area of 664 arcmin^, a factor of '--^ 1.14 larger 
than the total area surveyed by the sample used in the present work (~583 arcmin^). Their 90% completeness levels are in the 
range 22.0-22.4 mag at 3.6 /im and 4.5 /im, about a magnitude shallower than the 75% completeness levels. 

The stellar masses were derived by fitting the observed SEDs with a grid of models created with the PEGASE code (Fioc & 
Rocca-Volmerange 1997). The star-formation history (SFH) was parameterized with an exponentially declining SFH with an 
e-folding timescale r € [0.001, 100] Gyr, with allowed ages from 1 Myr to the age of the universe at the object's redshift. Seven 
discrete values of the metallicity were used, from Z = 0.005 Zq to Z = 5.0 Z©. A Salpeter (1955) IMF truncated at 0.1 and 100 
Mq was adopted. Dust extinction was modeled using the Calzetti et al. (2000) curve with Ay G [0,5]. While other SED-modeling 
assumptions were used to test how the stellar masses changed by changing IMF, extinction curve, stellar population synthesis 
model, and star formation history, the systematic effects of these changes on the derived SMFs were not explicitly quantified nor 
discussed in Perez-Gonzalez et al. (2008). 

SMFs were derived by integrating the bivariate luminosity-stellar mass function (the estimation of which was performed with 
a stepwise maximum likelihood technique) over all luminosities. The resulting SMFs were then fitted with a Schechter (1976) 
function. The low-mass end of the SMF at z > 1.6 was constrained by combining their results with other estimates of the SMFs 
from the Uterature. The redshift-dependent completeness limit were derived assuming to be complete for stellar masses larger 
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Fig. C5. — Comparison between the SMFs from Perez-Gonzalez et al. (2008) and from this work at 3.0 < z < 4.0 (left panels), 2.0 < z < 3.0 (middle panels), 
and 1.3 < z < 2.0 (right panels). Symbols as in Fig. CI. The SMFs of Perez-Gonzalez et al. (2008) have been scaled by -0.2 dex along the A:-axis to take into 
account difference in the adopted IMF. Note that the SMFs from Perez-Gonzalez et al. (2008) are actually derived at 3.0 < z < 3.5 and 3.5 < z < 4.0 (left panels; 
filled black circles and squares, respectively); 2.0 < z < 2.5 and 2.5 < z < 3.0 (middle panels; filled black circles and squares, respectively); 1.3 < z < 1-6 and 
1.6 < z < 2.0 (right panels; filled black circles and squares, respectively). 

than the stellar mass corresponding to a passively evolving stellar population formed at z ~ cx) with no extinction and having 
a 3.6 /im flux equal to the 75% completeness level of the IRAC sample. Only galaxies with stellar masses larger than this 
completeness level were used in their analysis, and no completeness correction was carried out to recover the SMF at smaller 
masses. 

A direct comparison between the SMFs of Perez-Gonzalez et al. (2008) and the SMFs measured in our work is shown in 
Fig. C5. Note that the redshifts bins used in Perez-Gonzalez et al. (2008) are smaller than the redshift intervals used in our work, 
i.e., our redshift intervals are further split in two in Perez-Gonzalez et al. (2008). 

As shown in Fig. C5, the SMFs from Perez-Gonzalez et al. (2008) are generally in good agreement with the SMFs derived in 
our work, at all redshifts, especially at the low-mass end. Only at the very high-mass end in the redshift ranges 2.0 < z < 3.0 
and 1 .3 < z < 2.0, the number densities derived from Perez-Gonzalez et al. (2008) look slightly larger than those derived in our 
work. These differences are only barely significant, and are definitely not statistically significant at all once cosmic variance and 
systematic uncertainties are included in the total error budget. We therefore conclude that the SMFs derived from Perez-Gonzalez 
et al. (2008) are in good agreement with the SMFs derived in our work, with the latter better sampling the low-mass end of the 
SMFs by 0.3-0.6 dex in stellar mass. 



REFERENCES 

Allen, C. W. 1976, Astrophysical Quantities. University of London, The Avni, Y, & Bahcall, J. N. 1980, ApJ, 235, 694 

Athlone Press, 264 Baugh, C. M., 2006, RPPh, 69, 3101 

Amouts, S., Cristiani, S.,Moscardini, L.,Matarrese, S., Lucchin, F.,Fontana, Bell, E. F., Mcintosh, D. H., Katz, N., & Weinberg, M. D. 2003, ApJS, 

A., Giallongo, E. 1999, MNRAS, 310, 540 133149, 289 



The SMFs at 1 .3 < z < 4.0 and their uncertainties 



33 



Benitez, N. 2000, ApJ, 536, 571 

Blain, A. W., Jameson, A., Small, I., Longair, M. S., Knelb, J.-P., & Ivlson, 

R. J. 1999, MNRAS, 309, 715 
Blain, A. W., Jameson, A., Small, I., Longair, M. S., Knelb, J.-P., & Ivlson, 

R. J. 1999, MNRAS, 309, 715 
Blain, A. W., SmaU, I., Ivlson, R. J., & Knelb, J.-P. 1999, MNRAS, 302, 632 
Blanton, M. R., & Rowels, S. 2007, AJ, 133, 734 
Bolzonella, M., Miralles, J.-M., Pell6, R. 2000, A&A, 363, 476 
Borch, A., et al. 2006, A&A, 453, 869 

Bouchet, P., Lequeux, J., Maurice, E., Prevot-Bumiclion, M. L. 1985, A&A, 
149, 330 

Brammer, G. B., van Dokkum, R G., & Coppi, P. 2008.ApJ, 686, 1503 

Bruzual, G., & Chariot, S. 2003, MNRAS, 344, 1000 

Bruzual, G. 2007, in Vazdekls A., Peletier R. F., eds, lAU Symp. Vol. 241, 

On TP-AGB Stars and the Mass of Galaxies. Cambridge Univ. Press, 

Cambridge, p. 125 
Bundy, K., et al. 2006, ApJ, 651, 120 

Cassisi, S., Castellani, M., & Castellani, V. 1997, A&A, 317, 108 

Chabrier, G. 2003, PASP, 115, 763 

Chariot, S., & Bruzual, G. 2008, in preparation 

Calzetti, D., Armus, L., Bohlin, R. C, Kinney, A. L., Koomneef, J., & 

Storchi-Bergmann, T. 2000, ApJ, 533, 682 
Cole, S., et al. 2001, MNRAS, 326, 255 

Coleman, G. D., Wu, C.-C, Weedman, D. W. 1980, ApJS, 43, 393 
Conselice, C. J., Blackbume, J. A., & Papovlch, C. 2005, ApJ, 620, 564 

Daddi, E., et al. 2007, ApJ, 670, 173 
Dave, R. 2008, MNRAS, 385, 147 

De Lucia, G., Kauffmann, G., & White, S. D. M. 2004, MNRAS, 349, 1 101 
De Lucia, G., & Blaizot, J. 2007, MNRAS, 375, 2 

Dickinson, M., Papovich, C, Ferguson, H. C, & Budavari, T. 2003, ApJ, 
587, 25 

Drory, N., Bender, R., Feulner, G., Hopp, U., Maraston, C, Snigula, J., & 

Hill, G. J. 2004, ApJ, 608, 742 
Drory, N., Salvato, M., Gabasch, A., Bender, R., Hopp, U., Feulner, G., & 

Pannella, M. 2005, ApJ, 619, LI 11 
Efstathiou, G., Ellis R. S., & Peterson, B. A. 1988, MNRAS, 232, 431 
Eisner, F., Feuhier, G., & Hopp, U. 2008, A&A, 477, 503 
Fagotto, F, Bressan, A., Bertelll, G., & Chiosl, C. 1994, A&AS, 104, 365 
Fazio, G. G., et al. 2004, ApJS, 154, 10 
Fioc, M., & Rocca-Volmerange, B. 1997, A&A, 326, 950 
Fontana, A., et al. 2003, A&A, 594, L9 
Fontana, A., et al. 2004, A&A, 424, 23 
Fontana, A., et al. 2006, A&A, 459, 745 

Fontanot, P., De Lucia, G., Monaco, P, SomervlUe, R. S., Santlni, P. 2009, 

MNRASaccepted [arXiv:0901.1130] 
Forster Schreiber, N. M., et al. 2006, AJ, 131, 1891 
Franx, M., et al. 2003, ApJ, 587, L79 

Franx, M., van Dokkum, P. G., Forster Schreiber, N. M., Wuyts, S., Labb6, 

1., Toft, S. 2008, ApJ, 688, 770 
Gawiser, E., et al. 2006, ApJS, 162, 1 
Gehrels, N. 1986, ApJ, 303, 336 

Giallongo, E., Salimbeni, S., Menci, N., Zamorani, G., Fontana, A., 

Dickinson, M., Cristiani, S., Pozzetti, L. 2005, ApJ, 622, 116 
Giavalisco, M., et al. 2004, ApJ, 600, L93 
Glazebrook, K., et al. 2004, Nature, 430, 181 
Grazian, A., et al. 2006, A&A, 449, 951 
Hauschildt, P H., Allard, P., & Baron, E. 1999, ApJ, 512, 377 
Heidt, J., et al. 2003, A&A, 398, 49 
Hora, J. L., et al. 2004, Proc. SPIE, 5487, 77 
llbert, O., et al. 2006, A&A, 457, 841 
lovino. A., et al. 2005, A&A, 442, 423 
Kauffmann, G., & Heahnelt, M. 200, MNRAS, 311, 576 



Kendall, M. G., & Stuart, A. 1961, The Advanced Theory of Statistics, Vol. 2, 

Griffin & Griffin, London 
Kinney, A. L., Calzetti, D., Bohlin, R. C, McQuade, K., Storchi-Bergmann, 

T., & Schmitt, H. R. 1996, ApJ, 467, 38 
Kriek, M., et al. 2007, ApJ, 669, 776 
Kriek, M., et al. 2008, ApJ, 677, 219 
Kroupa, R 2001, MNRAS, 322, 231 
Labb6, 1., et al. 2003, AJ, 125, 1107 
Larson, R. B. 2005, MNRAS, 359, 211 
Lawrence, A., et al. 2007, MNRAS, 379, 1599 
Le Pevre, O., et al. 2005, A&A, 439, 845 
Maraston, C. 2005, MNRAS, 362, 799 

Maraston, C, Daddi, E., Renzini, A., Cimatti, A., Dickinson, M., Papovich, 
C, Pasquali, A., & Pirzkal, N. 2006, ApJ, 652, 85 (Erratum: 2007, ApJ, 
656, 1241) 

Marchesini, D., et al. 2007, ApJ, 656, 42 

Marchesinl, D., & van Dokkum, P. 2007, ApJ, 663, L89 

Marigo, P, & Girardi, L. 2007, A&A, 469, 239 

Monaco, P, Pontanot, P, & Taffoni, G. 2007, MNRAS, 375, 1189 

Muzzin, A., et al. 2009, ApJ, submitted 

Papovich, C, Dickinson, M., Ferguson, H. C. 2001, ApJ, 559, 620 
Papovich, C, et al. 2006, ApJ, 640, 92 
Perez-Gonzalez, P G., et al. 2008, ApJ, 675, 261 
Pozzetti, L., et al. 2007, A&A, 474, 443 

Prevot, M. L., Lequeux, J., Prdvot, L., Maurice, E., Rocca-Volmerange, B. 

1984, A&A, 132, 389 
Quadri, R., et al., 2007, AJ, 134, 1 103 
Reach, W. T., et al. 2005, PASP, 1 17, 978 
Reddy, N. A., & Steidel, C. C. 2009, ApJ, 692, 778 
Rudnick, G., et al. 2003, ApJ, 599, 847 
Rudnick, G., et al. 2006, ApJ, 650, 624 
Salpeter, E. E. 1955, ApL 121, 161 

Sandage, A., Tammann, G. A., & Yahil, A. 1979, ApJ, 232, 352 
Scarlata, C, et al. 2007, ApJS, 172, 494 
Schechter, P 1976, ApJ, 203, 297 

Schlegel, D. J., Pinkbeiner, D. P, & Davis, M. 1998, ApJ, 500, 525 
Schmidt, M. 1968, ApJ, 151, 393 

SomervlUe, R. S., Hopkins, P. P., Cox, T. J., Robertson, B. E., Hemquist, L. 

2008, MNRAS, 391, 481 
SomervlUe, R. S., & Primack, J. R. 1999, MNRAS, 310, 1087 
SomervlUe, R. S., Primack, J. R., & Faber, S. M. 2001, MNRAS, 320, 504 
Spergel, D. N., et al. 2007, ApJS, 170, 377 

Springel, v.. White, S. D. M., Tormen, G., & Kauffmann, G. 2001, MNRAS, 
328, 726 

Springel, v., et al. 2005, Nature, 435, 629 
Taylor, E. N., et al. 2008, ApJ, submitted 
Toft, S., et al. 2007, ApJ, 671, 285 
van Dokkum, R G., et al. 2006, ApJ, 638, 59 

van Dokkum, P G. 2008, ApJ, 674, 29 

van Dokkum, P G., et al. 2009, PASP, 121, 2 

Vergani, D., et al. 2008, A&A, 487, 89 

Wang, J., De Lucia, G., Kitzbichler, M. G., White, S. D. M. 2008, MNRAS, 
384, 1301 

WiUdns, S. M., Hopkins, A. M., Trentham, N., Tojeiro, R. 2008, MNRAS, in 

press 

Wuyts, S., et al. 2007, ApJ, 655, 51 

Wuyts, S., Labbee, 1., Forster Schreiber, N. M., Franx, M., Rudnick, G., 

Brammer, G. B., & van Dokkum, P G. 2008, ApJ, 682, 985 
Zucca, E., et al. 2006, A&A, 455, 879 



